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SECTION  AA 


Figure  1.  Schpm.’Eic  of  the  Mixing  Layer 


Nav  ier-Stokes  equations  and  the  two-equation  turbulence  model 
(for  the  turbulence  kinetic  energy  k  and  its  rate  of  dissipation 
e  )  ,  the  numerical  simulation  of  the  Stanford  mixing  layer  became 
of  interest. 

Our  computational  investigations  have  been  examining  tur¬ 
bulent  recirculating  flowfields  characteristic  of  the  CBCC  (see, 
e.g.,  Kr  ishnamurthy ,  et  al.°)  by  means  of  the  TEACH1-0  code  which 
involves  elliptic  flowfields.  The  mixing-layer  development, 
however,  is  an  example  of  an  evolving  flow  (without  flow  recir¬ 
culation)  which  is  governed  by  parabolic  equations.  Strictly 
speaking,  a  numerical  modeling  based  upon  an  elliptic  formulation 
is  nit  necessary  to  address  the  mixing-layer  problem.  Although 
the  experimental  configuration19  does  not  involve  recirculating 
flow,  it  does  represent  a  confined  flowfield,  for  which  the  TEACH 
procedure  is  appropriate.  The  availability  of  comprehensive  data 
t  ro.n  the  Stanford  mixing  layer,  therefore,  constitutes  a  useful 
means  of  model  validation. 

(  2  )  Expe  r  iraen  t  a_l_  Details 

As  seen  schematically  in  Figure  1,  the  Stanford  mixing-layer 
experiment  involves  a  flow  facility  which  consists  of  a 
4-cm-high,  10-cm-wide  and  20-cm-lony  test  section.  Upstream  of 
the  test  section  is  a  we  1 1 -des ig ned  two-dimensional  contraction 
in  which  is  located  a  splitter  plate,  midway  between  the  con¬ 
fining  walls  at  the  top  and  bottom.  A  mixing  layer  is  formed 
downstream  of  the  trailing  edge  of  the  splitter  plate,  with  un- 
e  jue  l  velocities  in  the  initially  separate  upper  and  lower 
streams.  Nitrogen  is  the  carrier  gas  in  both  streams  flowing 
continuously  at  free  stream  velocities  of  up  to  20  m/ s .  Reacting 
flow  experiments  have  employed  dilute  concentrations  of  NO  in  one 
stream  and  O2/O3  in  the  other.  The  hot-wire  anemometry 
measurements  of  present  interest  were  made  in  nonreacting  flows 
in  the  x~z  plane  located  midway  along  the  width  of  the  test 
e-ction  ( Figure*  l).  The  available  lata  correspond  t  1  initial 
velocities  of  h  m/s  in  the  upper  and  3  m/s  in  the  1  >wer  streams. 


SECTION  II 

REYNOLDS- AVERAGED  COMPUTATIONS 


This  section  deals  with  the  numerical  modeling  of  isothermal 
turbulent  flowfields  through  the  framework  of  Reynolds-averaged 
Navier-Stokes  equations.  Included  are:  the  predictive  results 
of  the  two-dimensional  mixing  layer  and  a  comparison  of  these 
with  the  experimental  data;  the  numerical  modeling  of  the  CBCC 
flowfields;  and  a  description  of  the  predictions  of  the 
vortex-center  characteristics  in  the  CBCC. 

1.  CONFINED  TWO-DIMENSIONAL  MIXING  LAYER 

The  numerical  simulation  of  the  Stanford  mixing-layer 
experiment^  is  described  below. 

a  •  I  nt  r oduct  i <5 n 

The  plane  turbulent  mixing  layer  formed  downstream  of  a 
splitter  plate  separating  two  uniform  and  initially  parallel 
streams  represents  one  of  the  simplest  and  most  extensively 
studied  flows.  The  most  basic  flowfield  in  this  configuration  is 
the  classic  mixing  layer  arising  from  two  streams  of  unequal 
velocities  but  of  the  same  species.  Apart  from  its  intrinsic 
interest,  such  a  mixing  layer  is  often  adopted  in  experimental 
studies  of  turbulent  mixing  in  nonreacting  and  reacting  flows. 

The  reomt  absorption  measurement^  of  a  conserved  scalar  in  a 
re  act  in j  flow  under  fast  chemistry  (due  to  two  initially  unmixed 
reactants  of  trace  concentration  present  in  inert  carrier 
streams)  is  a  good  example. 

(  1 )  Bac  kground 

Extensive  hot-wire  anemometry  data  of  the  mean  and  fluc¬ 
tuating  velocity  components  in  nonreacting  flows  (when  one  or 
both  of  the  reactant  species  are  absent )  were  available  from  the 
Stanford  University  experiment. ^  In  view  of  our  >ngoing 
t  lowtield  computat  ions  employing  the  Reynolds-  a  v<.*  raged 
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behavior  of  the  CBCC.  However,  within  the  scope  of  the  present 
program,  that  investigation  was  not  made.  Instead,  attention  was 
focused  on  the  influence  of  boundary  conditions,  numerical  para¬ 
meters,  and  computational  grid  on  the  CBCC  flowfield  predictions 
when  only  the  annular  stream  is  present. 

4.  OUTLINE  OF  REPORT 

The  Reynolds-averaged  computations  of  the  CBCC  flowfields  are 
reported  in  Section  II.  Section  III  first  discusses  briefly  the 
preliminary  results  of  the  time-dependent  calculations  which  led 
to  the  shedding-1 ike  behavior  and  then  reports  the  results  of  the 
comprehensive  examination  of  the  unsteady-flow  modeling.  The 
conclusions  of  the  present  modeling  research  program  and  the 
recommendations  in  light  of  these  conclusions  are  outlined  in 
Section  IV. 
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Dynamics  Laboratory  ( AFWAL/F IMM )  vectorized,  compressible, 
time-dependent  code.  Since  the  ability  of  the  time-dependent 
calculations  to  provide  physically  acceptable  solutions  of  the 
CBCC  flowfield  depended  on  the  successful  implementation  of  the 
boundary  conditions  appropriate  to  the  modeling  of  confined, 
subsonic  tlowfields,  this  subject  formed  the  main  aspect  of  the 
research.  The  examination  of  the  boundary  conditions  and  of  the 
question  of  acoustic  interactions  from  the  exit-plane  boundary 
was  to  be  made  for  two  simpler  flow  situations  of  the  CBCC. 

These  corresponded  to  the  CBCC  with  only  the  annular  flow  and 
with  only  the  central  jet.  Successful  implementation  of  the 
necessary  boundary  conditions  to  establish  the  dynamic  behavior 
in  the  two  cases  was  to  he  followed  by  the  numerical  investi¬ 
gation  of  the  CBCC  with  both  annular  and  central  flows  present. 

The  actual  scope  of  the  unsteady- f low  modeling  research, 
however,  entailed  a  comprehensive  examination  of  the  influence  of 
the  inflow  and  outflow  boundary  conditions  on  the  existence  of 
the  bluff-body  vortex  shedding  or  the  , auk  thereof.  This  is  easy 
to  see.  In  the  AFWAL/POSF  experiments  on  the  CBCC,  there  is  at 
present  no  reason  to  suspect  that  the  initial  and  boundary  con¬ 
ditions  are  time  dependent.  The  interior  flowfield,  however, 
does  exhibit  certain  time-dependent  features.  Indeed  the  prelim¬ 
inary  computations  with  the  AFWA.L/FI MM  code  appeared  to  indicate 
that  the  solution  of  the  time-dependent  governing  equations  pro¬ 
duced  a  t ime-depondent  flowfield  in  the  CBCC  interior.  That  such 
a  situation  may  not  be  unique  became  apparent  when  under  certain 
boundary  conditions  the  time-dependent  equations  failed  to 
generate  a  time-dependent  solution  in  the  interior  and,  in  fact, 
the  temporal  behavior  was  damped  out  after  several  time  steps. 

Of  course,  such  a  situation  is  common  when  the  physical  problem 
possesses  a  steady  solution  and  often  the  time-dependent  for¬ 
mulation  is  employed  to  yield  the  steady  solution  asymptotically 
for  large  times.  It  became  important,  therefore,  to  ascertain 
whether  such  was  the  case  with  the  CBCC  flowfield.  >n  this  con¬ 
text,  a  computat  iona  i.  investigation  involving  tii  -dependent 
boundary  conditions  would  have  been  helpful  to  taucidate  the 


flowfield.  Other  potential  contributors  of  perhaps  minor 
significance  were  thought  to  be  the  large-scale  structures  in 
the  boundary  layers  of  the  confining  duct  and  the  centerbody 
surface,  and  the  duct  acoustics. 

Irrespective  of  the  source  of  the  unsteadiness,  its  presence 
in  itself  r  ised  important  questions  with  respect  to  the 
numerical  modeling.  For  example,  is  the  CBCC  flowfield 
inherently  unsteady?  If  so,  how  relevant  and  useful  is  the 
prediction  based  on  Reynolds-averaged  formulation?  To  the  extent 
that  the  latter  appeared  to  show  reasonable  agreement  with  the 
time-averaged  measurements,  how  would  these  compare  with  the 
results  from  the  time  averaging  of  the  time-dependent  solutions? 
The  unsteady  flow  modeling  with  the  aid  of  time-dependent 
Nav ie r-Stokes  equations  was  expected  to  provide  answers  to  at 
least  some  of  these  questions. 

3.  SCOPE  OF  PRESENT  WORK 

The  work  reported  here  deals  with  the  time-averaged  and 
time-dependent  calculations  of  the  isothermal  flowfields  in  the 
CBCC.  The  former  involved  the  AFWAL/PORT  version  of  the  TEACH 
program  and  the  modifications  implemented  during  this  research. 
The  investigations  examined  the  large-  and  small-scale  CBCC's  and 
compared  the  predictions  with  the  experimental  data.  In  addition 
to  the  centerline  (axial)  distributions  and  the  radial  distrib- 
tions  (at  various  axial  locations)  of  the  mean  and  fluctuating 
fields,  the  computations  addressed  the  characteristics  of  the 
vortex  centers.  Also,  the  Reynolds-averaged  calculations 
examined  the  two-dimensional,  mixing-layer  experiment  at 
Stanford  University.  While  this  configuration  is  not  directly 
relevant  to  the  CBCC  flowfields,  this  numerical  simulation  served 
to  test  several  modeling  aspects.  No  computational  investigation 
with  the  improved  TEACH  code  recently  available  to  the  Air 
Forced  was  made  within  the  scope  of  this  program. 


The  time-dependent  computations  involved  the  axisymmetric 
version  of  the  Air  Force  Wright  Aeronautical  Laboratories,  Flight 


be  adequately  described  and  understood  in  terms  of  a  time-mean 
flow."  The  implications  of  Richardson's  suggestion,  that  the 
predictions  must  be  based  on  an  unsteady-flow  model  for  com¬ 
parison  with  and  interpretation  of  experimental  data  in  the  near¬ 
wake  region,  were  discussed  earlier^  for  the  CBCC  flowfields. 

No  time-dependent  predictions  of  the  CBCC  were  available, 
however,  for  an  examination  of  the  validity  of  the  Reynolds- 
averaged  predictions  of  the  flowfield. 


2.  OBJECTIVES 

Accordingly,  the  numerical  modeling  of  the  centerbody 
combustor  flowfields  was  considered  with  twofold  objectives. 
First,  the  Reynolds-averaged  computations  with  the  TEACH  code 
addressed  additional  refinements  in  CBCC  predictions  with 
improved  physical  and  numerical  models.  Improvements  in  physics 
were  concerned  with  the  modifications  to  the  standard  k-e  model 
to  account  for  streamline  curvature  and  the  preferential 
influence  of  normal  stresses  in  the  dissipation  equation.  The 
refinements  in  the  area  of  numerics  dealt  with  the  differencing 
schemes  and  included  the  incorporation  of  power-law  differencing 
schemes  and  modifications  of  the  hybrid-upwind  differencing 
schemes . 

While  the  steady-state  predictions  were  expected  to  show 
better  agreement  with  the  measured  results  because  of  the 
improvements  in  the  physical  and  numerical  models,  a  crucial 
aspect  missing  in  the  Reynolds-averaged  computations  was 
accounting  for  the  observed  dynamic  behavior  of  the  CBCC 
flowfield.  Thus,  the  second  objective  of  this  research  was 
concerned  with  the  computational  investigation  of  the  unsteady 
flows.  It  seemed  that  the  flowfield  unsteadiness  could  arise 
from  several  possible  causes.  Thus,  the  vortex  shedding  off  the 
centerbody  trailing  face,  large-scale  entrainment  structures  in 
the  shear  layer  between  the  annular  and  central  streams,  and  the 
coherent  structures  in  the  central  jet  could  all  ither  singly  or 
in  combination  impart  temporal  characteristics  to  the  CBCC 
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The  Reynolds-averaged  predictions  of  the  CBCC  flowfields 
appeared  to  provide  reasonable  approximations  to  the  experimental 
results,  especially  for  the  mean  flowfields  under  isothermal  con¬ 
ditions.  The  quantitative  accuracy  of  the  predictions  of  the 
fluctuating  fields  was  generally  poor.  Several  possible  sources 
for  this  discrepancy  have  been  identified,  as  discussed  below. 

The  TEACH  procedure  utilizes  the  hybrid-upwind  differencing 
scheme  for  the  convective  terms. 3 '12  This  scheme  is  known  to 
suffer  from  numerical  diffusion  problems. 6  The  CBCC  flowfield  is 
characterized  by  large  streamline  curvature  and  is  a  good  example 
of  the  complex  turbulent  flow  defined  by  Bradshaw, 13  due  to  the 
extra  strain  rates  associated  with  the  streamline  curvature.  It 
is  known  that  the  turbulent  shear  stress  and  the  degree  of  an¬ 
isotropy  between  the  normal  stresses  are  very  sensitive  to  cur¬ 
vature.  Since  the  standard  "k-e  model"  does  not  account  for 
streamline  curvature  effects,  this  was  considered  to  be  partly 
responsible  for  the  discrepancy  between  the  predicted  and  the 
measured  recirculation  lengths.  Moreover,  as  noted  by  Leschziner 
and  Rodi,14  the  dissipation  (e)  equation  in  the  standard  k-e 
model  requires  modification  in  recirculating  flows  to  account 
for  the  preferential  influence  of  the  normal  stresses  on  the  tur¬ 
bulence  transport.  Finally,  the  CBCC  flowfield  appeared  to 
exhibit  certain  flow  unsteadiness  associated  with  the  bluff- 
body  near  wake. 

The  applicability  of  Reynolds-averaged  description  to 
flowfields  possessing  large-scale  unsteadiness  has  been  a  subject 
of  much  debate.  It  is  known  that  the  separated  shear  layers  and 
the  recirculation  regions  downstream  of  a  bluff  body  are 
inherently  unsteady.  Indeed,  high-speed  photographs  of  com¬ 
bustion  in  the  CBCC^  have  revealed  the  propagation  of  large- 
scale  toroidal  vortices,  with  successive  flame  bursts  and 
relative  quiescence  of  random  nature.  As  pointed  out  by 
Richardson,  in  a  discussion  of  the  results  of  the  axisymmetric 
wake  behind  a  disk  by  Carmody,16  "it  appears  unlik.ly  that  a  flow 
possessing  distinct  and  large-scale  periodic  characteristics  can 
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SECTION  I 


INTRODUCTION 

This  final  report  documents  the  computational  fluid  dynamic 
research  performed  for  the  Air  Force  Wright  Aeronautical 
Laboratories,  Aero  Propulsion  Laboratory  (AFWAL/PO) ,  by  the 
University  of  Dayton.  The  computational  investigations  had  two 
overall  objectives:  (a)  conducting  the  steady-state  predictive 
modeling  of  the  confined  turbulent  mixing  in  the  near-wake  region 
of  the  AFWAL/POSF  centerbody  combustor  configuration  and  (b) 
performing  the  time-dependent  computations  of  the  centerbody 
combustor  flowfields. 

1 .  BACKGROUND 

The  confined  recirculating  flow  involving  dual  coaxial 
streams  in  the  near-wake  region  of  an  axisymmetric  bluff  body  has 
been  the  focus  of  an  ongoing  diagnostic  and  predictive  research 
program  at  the  Aero  Propulsion  Laboratory  (APL) .  A  1/5-scale 

model  of  the  APL  centerbody  combustor  configuration  (CBCC)  is  in 
operation  at  the  University  of  California,  Irvine  (UCI).4 

Numerical  predictions  of  the  CBCC  flowfields  have  been  the 
subject  of  several  research  studies  in  the  recent  past.3~9  These 
studies  dealt  with  the  steady-state  flowfield  computations  based 
on  the  Reynolds-averaged  Nav ier-Stokes  equations.  The  numerical 
calculations  employed  the  well  known  Teaching  Elliptic 
Axisymmetric  Characteristics  Heur ist ica lly  (TEACH) 10  procedure. 
The  modeling  of  turbulence  necessitated  by  the  Reynolds  averaging 
was  based  on  an  eddy  viscosity  computed  by  the  two-equation  model 
for  the  turbulence  kinetic  energy  k  and  the  rate  of  its  dissipa¬ 
tion  e .  1 1 
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For  these  conditions,  the  nonuniformities  in  the  inlet  mean 
velocities  were  less  than  1  percent.  The  streamwise  turbulence 
intensities  (u/U)  in  the  inlet  were  below  0.5  percent.  The 
z-profiles  of  the  streamwise  components  of  the  mean  (U)  and  the 
rms  (u)  velocities  were  measured  at  the  four  streamwise  stations 
of  5,  12,  15,  and  18  cm  downstream  of  the  trailing  edge  of  the 
splitter  plate . 

( 3 )  Scope  of  Computations 

The  predictive  calculations  have  employed  the  standard 
features  of  the  TEACH-type  numerics,  such  as  the  "hybrid"  upwind 
d  i  f  f  erenc  ingl  2  and  the  k- e  turbulence  modelH  (which  uses  a 
constant  value  of  0.09  for  the  parameter  cp).  Additional 
features  considered  in  the  present  calculations  are  the  power-law 
differencing  scheme^'^O  (which  retains  the  diffusive  effects  for  a 
larger  cell-Peclet  number  range,  viz.,  -10  <  Pe  <  10 ,  than  the 
"hybrid"  upwind  scheme)  and  the  streamline-curvature  correction^ 
which  introduces  a  curvature-dependent  (and  thus  nonconstant) 

Cy  in  the  k- e  model.  The  corrections  for  the  effect  of 
streamline  curvature  are  obtained  with  and  without  a  modification 
of  the  G-equat ion. 14  Present  computations  also  examine  the 
influences  of  nonuniform  inlet  profiles  and  the  exit-plane  loca¬ 
tion  of  the  computational  domain.  The  sensitivity  of  the  numeri¬ 
cal  predictions  to  different  aspects  of  the  modeling  is  noted  by 
a  comparison  of  predicted  and  measured  results. 

b .  Numerical  Computations 

The  TEACH  code  describes  a  finite-difference  computational 
procedure  to  solve  the  Reynolds-averaged  Navier-Stokes  equations. 
This  procedure  involves  the  primitive  (pressure  and  velocity) 
variables  instead  of  the  stream  f unct ion-vort ic ity  formulation 
and  thus  entails  a  direct  solution  of  the  velocity  and  pressure 
fields.  A  special  procedure  called  the  SIMPLE21  algorithm  is 
employed  for  this  purpose.  Further  details  concerning  the 
underlying  theory  and  the  computational  procedure  -re  available 
in  References  3  and  10. 
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(  1  )  Gov ern  in t[  Ki[uat  i o n_s 


For  the  mixing-layer  problem  o£  interest  here,  the  governing 
equations  are  formulated  in  the  planar  two-dimensional  geometry. 

H  and  V  are  the  time-mean  velocity  components  in  the  longitudinal 
(x)  and  the  transverse  (z)  directions.  Present  calculations  deal 
with  the  isothermal,  nonreacting  flow  of  the  carrier  (nitrogen) 
gas.  Therefore,  it  is  not  necessary  to  consider  the  species  and 
energy  conservation  equations.  The  only  other  dependent 
variables  of  interest  are  k  and  e.  The  governing  equations  for 
all  the  dependent  variables  can  be  expressed  in  the  general  form 

h  (pu  ^  +  h  (pv  +>  ~  h  (r*  fl*  ‘  h  {%  =  (1) 

where  denotes  any  dependent  variable  (time-mean  value).  In 
Eq.(l)  S,},  is  the  source  term  for  the  variable  <J>  which  includes 
true  source  terms  (such  as  those  due  to  chemical  reactions)  as 
well  as  the  terms  not  covered  by  the  first  four  terms 
(representing  the  convective  and  diffusive  contributions).  is 

the  effective  exchange  coefficient  for  the  transport  of  the 
variable  and  is  given  by 


r  ij)  -  Pe  f  f  /  °  <j>  (2) 

where  ueff  is  the  effective  viscosity  in  the  flowfield  and  is 
the  appropriate  effective  Prandt  1/Schmidt  number  for  each  . 

Table  1  summarizes  the  relevant  information  for  all  the  dependent 
variables. 


The  effective  viscosity  ue££  appearing  in  Eq.  (2)  is  given 
by 


vjeff  =  u  +  ut 


(3) 


whore  n  is  the  Laminar  viscosity  and  |it  is  the  turbulent  eddy 
viscosity.  The  latter  is  obtained  from 
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where  c y  is  usually  taken  to  be  a  constant  equal  to  0.09.  The 
above  procedure  involves  the  introduction  of  two  partial  dif¬ 
ferential  equations  for  k  and  e  which  are  solved  together  with 
the  conservation  equations  for  mass,  momentum,  and  energy.  We 
note  that  the  inclusion  of  these  additional  equations  has  been 
anticipated  in  the  general  formulation  of  Eq.  (1)  and  in  Table  1 
The  "standard"  k-e  model11  does  not  account  for  streamline  cur¬ 
vature  effects.  Following  Leschziner  and  Rodi,1^  ad  hoc  modifi¬ 
cations  have  been  introduced  by  us  to  incorporate  curvature 
effects  in  the  flowfield  modeling  of  the  centerbody  con¬ 
figuration. 3  Although  streamline-curvature  effect  is  not 
expected  to  be  significant  in  the  mixing-layer  problem,  present 
calculations  examined  this  question  also. 

( 2 )  Computational  Details 

Figure  2  shows  the  computational  domain  and  the  grid-point 
distribution  initially  adopted  in  the  numerical  calculations. 

The  chosen  grid  consists  of  41  longitudinal  nodes  and  42  lateral 
nodes  with  grid  spacings  that  are  nonuniform  longitudinally  and 
uniform  laterally.  The  location  of  the  exit  boundary  at  22.5  cm 
from  the  trailing  edge  of  the  splitter  plate  is  arbitrary.  The 
experimental  details  of  the  region  downstream  of  the  20-cm-long 
test  section  are  not  known.  To  find  the  sensitivity  of  the 
results  to  the  location  of  the  exit  boundary,  calculations  were 
also  carried  out  with  a  different  grid  consisting  of  51  longitu¬ 
dinal  nodes  and  with  the  exit-plane  located  at  a  distance  of 
31.5  cm . 

In  the  TEACH  computational  procedure  all  the  salient 
locations  such  as  the  top  and  bottom  walls,  the  trailing  edge  of 
the  splitter  plate,  and  the  exit  boundary  are  located  midway 
between  the  adjacent  grid  nodes.  Moreover,  a  staggered-gr id 
arrangement  is  employed  in  the  formulation  of  the  finite- 
difference  equations.  All  the  dependent  variables,  except  U  and 
V,  are  referred  to  at  the  grid  nodes.  U  and  V  are  calculated  at 
locations  midway  between  the  grid  nodes. 


( 3 )  Boundary  Conditions 

Table  2  summarizes  the  boundary  conditions  employed  in  the 
numerical  calculations.  In  view  of  the  elliptic  formulation  of 
Eq.  (1),  the  boundary  conditions  are  prescribed  on  all  the 
boundaries  of  the  computational  domain.  The  top  and  bottom 
boundaries  are  the  confining  walls  and  the  dependent  variables  U, 
and  e  are  prescribed  by  wall-function  formulation.  V  vanishes  on 
these  boundaries.  The  right  boundary  is  the  outflow  boundary  at 
which  the  streamwise  derivatives  of  all  the  dependent  variables 
are  set  to  zero.  The  left  boundary  represents  the  inflow  boun¬ 
dary.  The  inlet  profiles  of  U  correspond  to  experimental  con¬ 
ditions  (e.g.,  Uin  is  6  ra/s  in  the  upper  stream  and  3  m/s  in  the 
lower  stream).  With  well  designed  inlets,  the  transverse 
velocity  V  is  close  to  zero  and  is  treated  as  such.  Numerical 
calculations  have  been  performed  for  both  uniform  and  nonuniform 
inlet  distributions  of  U. 

The  specification  of  k  and  e  at  the  inlets  requires  some 
elaboration.  The  inlet  profile  of  k  is  obtained  from 

k  =  TURBIN  x  Ujn  (5) 

where  U^n  denotes  the  mean  longitudinal  velocity  at  the  inlet  and 
TURBIN  is  the  parameter  (FORTRAN  variables  in  the  TEACH  code) 
which  remains  to  be  specified.  We  note  that  the  experimental 
data  on  the  turbulence  intensity  measurements  in  the  inlet 
provide  a  basis  for  selecting  the  appropriate  value  of  TURBIN. 
Furthermore,  in  the  absence  of  turbulence  intensity  measurements 
in  all  three  orthogonal  coordinate  directions,  the  determination 
of  k  from  the  intensity  results  in  one  or  two  directions  (e.g., 
the  Stanford  mixing-layer  data  provide  u/U  only)  requires  the 
assumption  of  isotropy. 

The  inlet  distribution  of  e  is  specified  through  the 
prescription  of  the  inlet  turbulence  length  scale,  according  to 
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where  i  is  the  inlet  turbulence  length  scale.  For  the  two- 
dimensional  mixing  layer,  l  is  taken  as 

l  =  Xh,  (7) 

where  X  is  a  parameter  (FORTRAN  variable  in  the  TEACH  code  is 
ALAMDA )  and  h  is  the  half  height  of  the  channel  (h  =  2  cm).  The 
value  of  0.56  for  X  used  in  the  present  calculations  is  based 
upon  earlier  studies.3'14  In  addition,  the  sensitivity  of  the 
predictions  has  been  checked  for  X  =  0.05. 

c .  Results  and  Discussion 

Present  numerical  investigations  have  considered  several 
aspects  of  the  mixing-layer  flowfield.  The  results  reported  here 
correspond  to  initial  velocities  of  6  and  3  m/s  in  the  upper  and 
Lower  streams  of  nitrogen.  The  Reynolds  number  Rex  [based  upon 
the  distance  x  from  the  trailing  edge  of  the  splitter  plate  and 
defined  as  Rex  =  0.5  (Ui  +  U2)  x/ v ,  where  v  is  the  kinematic 
viscosity]  considered  ranges  from  14,300  for  x  =  5  cm  to  52,000 
for  x  =  18  cm.  The  numerical  predictions  are  presented  here  and 
compared  with  experimental  data. 

(  1  )  Influence  of  Inlet  Velocity  Profile 

The  predicted  results  of  the  mean  and  rms  streamwise  velo¬ 
city  components  in  the  mixing  layer  are  shown  in  Figures  3 
through  6  for  the  four  x  locations  of  5,  12,  15,  and  18  cm 
respectively.  The  predictions  corresponding  to  both  uniform  and 
nonuniform  inlet  distributions  of  U  and  k  are  compared  with  the 
hot-wire  data.  Since  the  available  experimental  data  at  the 
inlet  sections  do  not  go  all  the  way  to  the  splitter  plate  and  to 
the  top  and  bottom  walls,  some  minor  arbi t ra r i ness  arises  con¬ 
cerning  the  rapid  changes  in  U  and  u  close  to  these  boundaries. 
While  the  rapid  decay  of  the  mean  velocity  to  zero  does  not  pose 
problems,  the  rapid  increase  first  and  then  the  rapid  decrease  of 
the  rms  velocity  to  zero  at  the  boundaries  cause  ome  difficulty 
in  prescribing  the  nonuniform  inlet  profile  tor  the  numerical 
oa  Leu l a  t ions . 
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An  examination  of  the  results  indicates  that  the  agreement 
between  the  predictions  and  measurements  is  much  better  tor  the 
mean  velocity  component  than  for  the  rms  component.  For  the  mean 
velocity  component,  the  predicted  results  with  uniform  inlet  pro¬ 
files  show  very  good  agreement  with  the  experimental  data. 
Calculations  with  nonuniform  inlet  profiles  underpredict  the 
measurements  at  the  four  streamwise  locations.  For  the  rms  velo¬ 
city  component,  the  predictions  show  that  the  overall  trend  in 
the  mixing  layer  conforms  to  the  experimental  trend. 

Quantitative  agreement  for  the  magnitude  and  location  of  the  peak 
turbulence  intensity,  however,  is  generally  poor.  We  note  that 
the  calculations  overpredict  the  turbulent  mixing  rate,  espe¬ 
cially  near  the  origin  of  the  mixing  layer.  Farther  downstream, 
predictions  tend  to  indicate  better  agreement  with  the  measured 
data  inside  the  mixing  layer  and  in  this  case  also  the  calcula¬ 
tions  with  uniform  inlet  velocity  profiles  are  clearly  superior 
to  those  employing  nonuniform  profiles.  The  poor  agreement  be¬ 
tween  the  measurements  and  predictions  outside  the  mixing  layer 
may  be,  in  part,  due  to  the  assumption  of  isotropy  in  the  numeri¬ 
cal  modeling  [u  is  computed  as  (2k/3)®*^],  the  validity  of  which 
is  not  established  in  the  experiments.  Also,  the  very  rapid  rise 
in  the  intensity  predictions  near  the  top  and  bottom  walls  is  a 
consequence  of  the  normalization  with  respect  to  the  local  mean 
velocity  U  (which  vanishes  at  the  walls). 

(  2  )  Jjn  f  L  uence  of  Turbulence  Length  _S c  al  e ,  Differencing 
Sch erne,  a n d_  St  ream  1_  in e  Curvat u re 

The  predictive  calculations  discussed  earlier  were  based  upon 
the  "hybrid"  upwind  differencing  scheme  and  the  "standard"  k-e 
model  without  the  correction  for  the  streamline  curvature.  The 
inlet  turbulence  length  scale  parameter  X  was  taken  as  0.56. 
Additional  computations  were  completed  to  investigate  the  effect 
of  the  changes  in  these  areas  of  modeling  on  the  predictions. 

The  calculations  showed  that  there  is  no  discernible  difference 
in  the  mean  velocity  profile  at  all  four  x  locations  between: 

(^)  X  =  0.56  and  X  =  0.05,  (b)  upwind  and  power-law  differencing 
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schemes,  and  (c)  the  results  with  and  without  the  correction  tor 
streamline  curvature. 

Figures  7  through  10  show  the  results  o£  the  turbulence 
intensity  profiles  at  the  four  streamwise  locations.  It  appears 
that  the  variations  in  both  the  inlet  turbulence  length  scale  and 
the  differencing  schemes  do  not  have  significant  influence  on  the 
rms  velocity  field  at  the  four  x  locations.  The  curvature 
correction  shows  no  discernible  influence  on  the  rms  velocity 
predictions  at  the  first  two  streamwise  locations.  Farther 
downstream  (at  x  =  15  and  18  cm),  however,  a  small  but  noticeable 
effect  is  seen  inside  the  mixing  layer. 

(  3  )  Cl  the  r  Effects 

All  the  foregoing  results  were  obtained  from  calculations 
that  employed  the  computational  domain  with  the  exit-plane 
located  at  x  =  22.5  cm.  Moreover,  the  streamline  curvature 
correction  did  not  involve  the  dissipation-equation  modification 
of  Reference  14.  Therefore,  a  parametric  examination  of  the 
effect  of  a  different  exit-plane  location  (x  =  31.5  cm)  and 
^-equation  modification  was  completed.  These  calculacions  showed 
n  >  difference  in  the  mean  and  rms  predictions  at  the  first  three 
x  locations.  Figure  11  which  presents  the  mean-velocity  profile 
at  x  -  18  cm  indicates  that  there  is  no  discernible  effect  inside 
the  mixing  layer.  The  turbulence  intensity  profile  at  x  =  18  cm 
seen  in  Figure  12  shows  a  very  small  effect  in  the  mixing  layer. 

d  .  C o nc  1  u sion s 

•  The  numerical  calculations  employing  the  TEACH-type  numerics 
and  the  k-c  turbulence  model  provide  physically  correct  predic¬ 
tions  of  the  two-dimensional,  isothermal  mixing  layer. 

•  The  numerical  predictions  of  the  mean  streamwise  velocity 
profile  at  different  downstream  locations  show  good  agreement 
with  the  hot-wire  data  from  the  Stanford  mi x i ng- layer  experiment. 


( 2 )  Turbu lence  Mode  1 

The  turbulent  eddy  viscosity  p^  (the  effective  viscosity 
Meff  appearing  in  the  governing  equations  is  given  by  ^  +  11 
where  p  is  the  laminar  viscosity)  in  the  k-e  model  is  obtained 
from  pt  -  Cypk^/e,  where  p  is  the  mass  density  and  cp  is  a 
constant  equal  to  0.09. 

Isotropy  is  assumed  in  obtaining  k  from  the  rms  velocity 
components,  according  to  k  =  (3/2)w2,  where  w  is  the  axial  rms 
velocity  component.  The  inlet  profile  of  k  is  given  by  0.03 
for  the  annular  jet  and  by  0.03  for  the  central  jet  (see  SS 
and  KP) .  The  required  inlet  profile  for  e  is  specified  as 
discussed  below. 

Inlet  Turbulence-Length  Scale.  As  in  KP ,  we  obtain  the 
inlet  profile  of  e  from  e  =  where  the  inlet  turbulence- 

length  scale  l  is  given  by  £  =  X5.  Here  X  is  a  specified 
constant  and  6  is  a  characteristic  reference  length  given  by 
( Dg-D) / 2  for  the  annular  jet  and  d/2  for  the  central  jet.  KP 
used  0.3333  and  0.5556  for  X.  The  former  value  is  equivalent  to 
0.03  used  in  3S  (since  they  employ  e  =  cvk^'^/H)  and  the  latter 
value  is  used  in  LR.  Present  results  are  based  upon  X  =  0.5556. 

Stream! ine  Curvature  Correction .  Curvature  modifications  to 
turbulence  models  have  been  attempted  by  a  number  of  researchers 
(e.g.,  see  LR;  Humphrey  and  Pourahmad i ; 3 1  Rodi  and  Scheuerer;32 
and  the  references  cited  in  them).  A  majority  of  these  modifica¬ 
tions  have  been  concerned  with  the  k-e  model.  Present  numerical 
investigation  has  introduced  a  curvature-dependent  (and  hence 
nonconstant)  c u  into  the  standard  k-e  model.  Following  LR,  we 
have 

w 
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c 

of  a  streamline,  s  and 
to  the  streamline,  and 


-  „  aw  W 

cu=  -K1K2/h+8Ki(k/c)  (-^r  +  R- 

c 

where  Rc  is  the  local  radius  of  curvature 
n  denote  the  coord inates  along  and  normal 
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TABLE  3 

FLOWFIELD  CONDITIONS 
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Figure  14.  Computational  Domain  for  CFSCC. 
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data2f30  to  elucidate  the  distinctions  between  the  small-  and 
large-scale  CBCC.  This  was  necessary  in  view  of  the  differing 
nature  of  the  predictions. 

(  2  )  Sc oj3 e  of  Computations 

This  study  deals  with  the  isothermal  flow-field  predictions. 
As  in  SS  and  KP ,  present  work  employs  the  TEACH  computational 
procedure  (Gosman  and  IderiahlO),  and  the  standard  features  of 
the  numerics  thereof.  An  additional  feature  considered  here  is 
the  power-law  differencing  scheme^  in  place  of  the  hybrid  upwind 
scheme. 32  while  both  the  APL  and  the  UCI  configurations  were 
studied,  only  a  limited  parametric  range  of  the  flowfield  con¬ 
ditions  investigated  is  reported. 

b .  N ume r ical  Computations 

The  details  concerning  the  application  of  the  TEACH  proce¬ 
dure  to  the  CBCC  flowfields  are  available  in  Kr i shnamurthy  et 
al.^  and  therefore  are  not  repeated  here. 

( 1 )  Computational  Details 

Figure  14  shows  the  present  computational  domain  and  the 
grid-point  distribution.  The  grid  consists  of  41  axial  nodes  and 
39  radial  nodes  with  a  nonuniform  spacing  in  both  directions  to 
ensure  adequate  spatial  resolution  in  flowfield  regions  with 
large  gradients  of  the  flow  variables.  The  calculations  with 
this  grid  require  about  170,0008  words  of  memory  on  the  CDC  CYBER 
computer.  The  converged  solutions  are  attained  in  less  than  900 
iterations,  with  each  iteration  taking  less  than  a  second  of  com¬ 
puter  time. 

The  results  herein  correspond  to  the  conditions  shown  in 
Table  3,  with  uniform  inlet  profiles  for  the  mean  axial  velocity, 
corresponding  to  W&  in  the  annular  jet  and  Wp  in  the  central  jet. 
The  Reynolds  number  Re^  of  the  air  flow  is  based  upon  D  and 
the  reference  velocity  Wd  in  the  duct.  The  Roynoils  number  Rep 
of  the  CO2  flow  is  based  upon  d  and  Wp. 
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SS  and  KP  employed  the  two-equation  model  for  the  turbulent 
kinetic  energy  k  and  its  dissipation  rate  e.H  This  model  in  its 
standard  version  does  not  account  for  streamline  curvature 
effects.  Its  predictions  overestimated  the  extent  of  the  recir¬ 
culation  region  in  the  APL  CBCC  compared  with  the  laser  Doppler 
anemometry  (LDA)  data  of  Lightman  et  al. 23-24  This  was  in 
contrast  to  the  experience  of  several  previous  studies  of  tur¬ 
bulent  recirculating  flows  where  significant  underprediction  was 
observed  (e.g.,  see  Pope  and  Whitelaw;25  Militzer  et  al.;26 
Gosman  et  al.;27  Durst  and  Rastogi28).  Also,  SS  and  KP  overpre¬ 
dicted  the  central  jet  CO2  mass  flow  rate  required  in  the  experi¬ 
ments  24  to  eliminate  the  centerline  reverse  air  flow.  Another 
aspect  in  which  the  prediction  was  deficient  concerned  the  axial 
distribution  of  the  centerline  rms  axial  velocity  component. 

It  is  known  that  the  turbulent  shear  stress  and  the  degree 
of  anisotropy  between  the  normal  stresses  are  very  sensitive  to 
streamline  curvature.29  Thus,  the  presence  of  large  curvature  in 
the  CBCC  may  have  contributed  in  part  to  the  discrepancy  between 
the  measurement  and  prediction.  Accordingly,  the  main  objective 
of  the  present  study  was  to  examine  the  influence  of  curvature 
corrections  in  the  prediction.  Recent  availability  of  extensive 
isothermal  LDA  data  from  UCI ,  as  well  as  LDA  (Lightman  et 
al.2)and  CO2  concentration  (Bradley  et  al.^O)  data  from  APL  has 
facilitated  this  inquiry. 

When  calculations  were  made  for  the  UCI  CBCC  with  the  stan¬ 
dard  k- e  model,  the  centerline  recirculation  length  was  underpre¬ 
dicted  in  line  with  the  experience  of  other  researchers. 

The  introduction  of  a  curvature  correction,  along  the  lines 
suggested  by  Leschziner  and  Rodi^  (hereafter  denoted  by  LR)  , 
resulted  in  only  a  partial  improvement  of  the  prediction. 
Therefore,  the  second  objective  was  to  examine  the  impact  of  the 
diffusion-equation  modification  considered  in  LR.  A  final  objec¬ 
tive  was  the  comparison  of  the  predictions  with  the  newer  APL 
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Figure  13.  Schematic  of  CBCC. 


32 


the  predicted  and  measured  results  of  the  axial  and  radial 
distributions  of  the  mean  and  rms  axial  velocity,  the  centerline 
stagnation  points,  and  the  axial  decay  of  the  centerline  CO2 
concentrat ion. 

a .  Centerbody  Combustor  Configuration 

The  centerbody  combustor  represents  a  complex  turbulent 
flowfield.  As  noted  by  Bradshaw, 13  a  complex  turbulent  flow  is 
one  which  cannot  be  predicted  with  acceptable  accuracy  by  methods 
developed  in  classical  thin  shear  layers.  The  assumptions 
inherent  in  the  thin-shear-layer  approximation  are  often 
invalidated  in  realistic  flowfields  by  the  presence  of  several 
features  (e.g.,  interacting  shear  layers  or  high  rates  of  strain 
associated  with  large  streamline  curvature).  Figure  13  shows 
schematically  the  large-  and  small-scale  CBCC.  This  involves 
confined  turbulent  mixing  of  an  annular  air  stream  and  a  central 
gas  (CO2  in  isothermal  experiments  and  C3H8  in  combusting 
experiments)  jet  in  the  near-wake  region  downstream  of  a  bluff 
body.  Figure  13  shows  that  the  ratio  D/d  is  rather  large  (29  for 
the  APL  and  23  for  the  UCI  configurations).  Indeed,  the  present 
interjet  separation  is  much  larger  than  that  studied  in  typical 
coaxial  jet  mixing  of  both  confined  and  unconfined  flowfields  in 
the  literature.  The  wide  separation  between  the  jets  and  the 
concomitant  presence  of  the  toroidal  recirculating  bluff-body 
wake  in  the  mixing  region  have  raised  some  interesting  flowfield 
ramifications  addressed  only  recently  in  numerical  predictions 
(e.g.,  see  Kr ishnamurthy ; 5  Sturgess  and  Syed,6  hereafter  denoted 
by  SS;  and  Kr ishnamurthy  et  al.,3  hereafter  denoted  by  KP) . 
Depending  on  the  strength  of  the  annular  and  central  jets,  the 
CBCC  flowfield  exhibits  wake-like  and  jet-like  characteristics 
under  isothermal  conditions  (e.g.,  see  Kr i shnamurthy ,  Wahrer,  and 
Cochran‘S).  Nevertheless,  the  juxtaposition  of  the  two  streams 
with  the  near-wake  region  in  a  confined  configuration  renders  the 
CBCC  flowfield  a  highly  "complex  turbulent  flow,"  in  the  sense  of 
Bradshaw.  ^  in  addition  to  the  interacting  she~r-  layers,  large 
streamline  curvature  is  present. 
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•  The  numerical  predictions  of  the  turbulence  (streamwise) 
intensity  are  consistent  with  the  experimental  trends  and  show 
fair  to  poor  quantitative  agreement. 

•  Both  the  mean  and  rms  velocity  predictions  are  not  signifi¬ 
cantly  affected  by  parametric  changes  in  the  numerical  modeling 
involving  inlet  turbulence  length  scale,  differencing  scheme, 
streamline  curvature  correction  in  the  k-e  model,  and  the  loca¬ 
tion  of  the  computational  domain  exit  plane. 

•  Predictions  based  on  uniform  inlet-velocity  profiles  show 
better  agreement  with  the  measurements  than  the  predictions  based 
on  nonuniform  profiles. 

•  The  mixing  layer  appears  to  exhibit  some  interactions  with 
the  boundary  layers  at  the  top  and  bottom  walls  (being  greater 
for  the  latter  than  for  the  former)  at  large  distances 
downstream.  These  interactions  may  explain,  in  part,  the  noted 
differences  between  the  predicted  and  measured  rms  velocity 
results,  when  there  is  excellent  agreement  for  the  mean  velocity. 
Indeed,  this  conclusion  is  consistent  with  the  recent 
observations  of  Wood  and  Bradshaw.22  They  report  significant 
changes  in  the  turbulence  structure  (with  the  mean-velocity 
profile  remaining  unaltered)  when  the  mixing  layer  is  influenced 
by  a  solid  surface.  while  their  mixing  layer  had  to  contend  with 
only  a  bottom  boundary  layer,  they  emphasize  that  "...similar, 
and  perhaps  even  larger,  changes  will  occur  in  a  two-stream 
mixing  layer  confined  by  a  tunnel  roof  as  well  as  the  floor." 

2.  LARGE-  ANI)  SMALL-SCALE  CBCC’S 

The  confined  turbulent  recirculating  flowfield  due  to  the 
isothermal  mixing  of  dual  coaxial  streams  in  the  near-wake  region 
of  the  CBCC  is  described  below.  The  finite-difference 
computations  examine  modifications  to  the  turbulence  model  to 
account  for  the  effects  of  streamline  curvature  and  for  the 
preferential  influence  of  normal  stresses  on  turbulence 
dissipation.  This  examination  is  facilitated  by  a  comparison  of 
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Influence  of  Curvature  and  Exit  Plane  Location  on 
the  rms  Velocity  Predictions  at  x  =  18  cm. 


t 


Figure  11 


Influence  of  Curvature  and  Exit  Plane  Location  on  the 
Mean  Velocity  Predictions  at  x  =  18  cm. 
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Figure  10.  Influence  of  Turbulence  Length  Scale,  Differencing 
Scheme  and  Curvature  on  the  rms  Velocity  Prediction 
at  x  =  18  cm. 
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Influence  of  Turbulence  Length  Scale,  Differencing 
Scheme  and  Curvature  on  the  rms  Velocity  Predictions 
at  x  =  15  cm. 
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Figure  8.  Influence  of  Turbulence  Length  Scale,  Differencing 

Scheme  and  Curvature  on  the  rms  Velocity  Predictions 
at  x  =  12  cm. 
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Influence  of  Turbulence  Length  Scale ,  Dif ferencincr 
Scheme  and  Curvature  on  the  rms  Velocity  Predictions 
at  x  =  5  cm. 


Figure  7 


K i  and  K2  are  constants  equal  to  0.267  and  -0.489  respectively. 
Although  an  element  of  arbitrariness  in  LR  was  noted  in  Reference 
3,  the  discussion  therein  emphasizes  the  utility  of  the 
formulation  of  LR  which  is  retained  here.  Also,  as  in  LR,  the 
calculations  impose  an  arbitrary  positive  lower  bound  of  0.026  on 


e  Modification.  We  employ  the  preferential  modification  of 
dissipation  proposed  in  LR  for  recirculating  flowfields. 

Although  additional  work  is  necessary  in  view  of  the  ad  hoc 
nature  of  this  modification,  its  inclusion  throws  some  light  on 
the  observed  distinction  between  the  small-  and  large-scale  CBCC. 

c .  Results  and  Discussion 

( 1 )  Pertinent  Results  from  SS  and  KP 

To  set  the  present  results  in  context,  we  show  in  Figure  16 
the  earlier  results  from  SS  and  KP.  With  identical  values  of  X 
(0.03  in  SS),  the  two  (grid  A  in  KP)  show  good  agreement  in  their 
degrees  of  underprediction  of  the  forward  stagnation  points  (FSP) 
and  overprediction  of  the  rear  stagnation  points  (RSP)  (the  FSP 
occurs  where  the  central  jet  loses  5 ' s  forward  momentum  entirely 
and  the  RSP  is  the  usual  end  of  the  recirculation  zone). 

The  measurements23-24  showed  that  for  W^  =  47  m/s  and  Wp 
=  135  m/s  (corresponding  respectively  to  an  air  mass  flow  of 
2  kg/s  and  CO2  mass  flow  of  16  kg/hr  in  the  APL  CBCC),  W  was  non¬ 
negative  on  the  centerline.  Indeed,  the  minimum  CO2  mass  flow 
which  eliminated  the  center 1 ine- f low  reversal  was  inferred  from 
the  measurements  to  be  14.7  kg/hr  in  SS.  The  predictions, 
however,  indicate  the  occurrence  of  the  centerline  reverse  flow 
(with  both  stagnation  points  present)  at  16  kg/hr. 

Curvature  Correction.  The  preliminary  results  in  KP  for 
streamline  curvature  correction  are  shown  for  two  CO2  flow  rates 
(viz.,  4  and  8  kg/hr).  It  must  be  noted  that  the  increase  in  X 
from  0.3333  to  0.5556  has  masked  the  effect  of  th  correction. 
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Figure  15.  Measured  and  Predicted  Centerline  Stagnation  Points 

for  APL  C3CC .  •,  ■  LDA  Results:  Lightman  et  al . 23,24 

The  Predictions  are  Based  upon  Hybrid  Upwind  Differ¬ 
encing  Scheme.  -  :  SS  (39  x  39),  X  =  0.3333;  O  : 

KP,  Grid  A  (41,  34),  X  =  0.3333; O:  KP ,  Grid  B 
(41,  34),  X  «  0.3333;  A:  KP ,  Grid  B  (41,  34), 

X  =  0.5556,  with  Curvature  Corrections.  Grid  B  had 
More  Nodes  in  the  Annular  Region  than  Grid  A. 


SS  and  KP  have  observed,  however,  that  the  increase  in  X  (without 
the  correction)  moves  monoton ica 1 ly  both  the  stagnation  points 
upstream,  thereby  resulting  in  greater  underprediction  of  the 
FSP.  With  the  correction,  the  FSP  moves  farther  downstream 
towards  the  measured  value.  The  degree  of  overprediction  of  the 
RSP  is  seen  greatly  reduced  with  the  larger  X  and  the  curvature 
correction.  Finally,  the  prediction  with  the  correction  has  elim 
inated  the  centerline-flow  reversal  at  the  CO2  flow  of  16  kg/hr. 

Unlike  the  good  agreement  for  the  centerline  stagnation 
points,  the  centerline  peak  negative  mean  axial  velocity  Wm  was 
overpredicted  with  the  curvature  correction.  Since  KP  had  not 
isolated  the  influence  of  X,  hybrid  upwind  differencing  (and  the 
propensity  for  numerical  diffusion  therein),  and  the  arbitrarily 
modified  grid  B,  it  was  essential  to  examine  the  effect  of  the 
curvature  correction  more  systematically. 

( 2 )  UCI  Configuration 

Figure  16  compares  the  present  predictions  with  the 
measurements.  We  note  that  the  standard  model  underpredicts  both 
the  stagnation  points  by  about  30  percent,  a  result  in  agreement 
with  other  studies.  Considerable  improvement  is  seen  when  the 
curvature  correction  to  c ^  is  introduced.  However,  the 
stagnation  points  and  Wra  are  still  underpredicted.  With  the 
e-equation  modification  in  addition  to  the  curvature  correction, 
excellent  agreement  is  seen  for  0  <  z/D  <  1.6.  Farther 
downstream  the  measured  recovery  is  greater  than  given  by  all 
the  predictions. 

Comparison  of  rms  Velocity.  The  results  seen  in  Figure  17 
do  not  confirm  the  good  agreement  observed  for  the  mean  axial 
velocity.  In  fact,  the  predictions  are  very  poor  especially  for 
z/D  <  1.3.  The  somewhat  better  agreement  seen  for  2  <  z/D  <  3  is 
not  consistent  with  the  greater  disparity  noted  for  the  mean 
velocity.  Of  course,  a  major  problem  in  the  predictions  is  the 
assumption  of  isotropy,  which  is  known  to  be  invalid  in  the  near 
wake.  Furthermore,  good  agreement  for  W  and  poor  agreement  for  w 
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Figure  17. 
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Measured  (Brum  et  al .  )  and  Predicted  Centerline  rms 
Axial  Velocity  for  UCI  CBCC . 


see 


appear  to  be  characteristic  of  complex  turbulent  flows  (e.g.. 

Wood  and  Bradshaw^ 

Radial  Distributions.  The  preiicted  radial  profiles  at  four 
axial  locations  in  Figures  18  and  19  include  both  the  c^  and  e 
modifications.  Good  agreement  is  seen  and  the  comparison  is  par¬ 
ticularly  striking  at  the  axial  location  well  within  the  recir¬ 
culation  region  (z/D  ~  0.17)  and  at  the  one  near  the  RSP 
( z/D  ~  1.31). 

Flowfield  Without  the  Central  Jet.  There  is  a  marked  dif¬ 
ference  in  Figure  20  where  the  central  jet  is  absent.  The  RSP  is 
underpredicted  by  about  10%  but  the  predicted  magnitude  and  the 
axial  location  of  wm  differ  by  a  factor  of  2  from  the  measured 
values.  The  rm.  component,  however,  shows  better  agreement  be¬ 
tween  the  prediction  and  the  measurement.  This  behavior  is 
clearly  different  from  that  of  the  nonzero  central  jet  seen  in 
Figures  16  and  17.  A  parametric  dependence  emerging  from  these 
observations  appears  to  emphasize  the  ad  hoc  nature  of  the  modi¬ 
fications  examined. 

(  3  )  API,  Configuration 

The  centerline  profiles  of  W  and  w  for  the  large-scale  CBCC 
are  seen  in  Figures  21  through  23  for  three  different  values  of 
Wp .  The  predictions  of  the  standard  k-e  model  and  those  with  the 
correction  for  curvature  (with  and  without  the  e  modification) 
are  compared  with  the  LDA  data  of  Lightman  et  al . ^ 

Measurements .  The  experimental  data  for  the  mean  velocity 
differ  from  the  earlier  results  (Lightman  et  al . ^ ^ ) . The  loca¬ 
tions  of  the  RSP  are  at  a  z/D  of  1  instead  of  0.9  (see  Figure 
15).  Figure  23  indicates  a  small  region  of  centerline-flow 
reversal  (0.75  <  z/D  s  0.975).  The  newer  data  thus  show  better 
agreement  with  the  predictions  in  SS  and  KP.  Indeed,  for  the 
case  of  zero  CO?  flow,  the  RSP  result  of  SS  even  shows  a  small 
underpreu  ict.  i  on  . 
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Figure  19.  Measured  (Brum  et  al.  )  and  Predicted  rms  Axial 
Velocity  Radial  Profiles  for  UCI  '"IBCC . 


45 


vM/« 

—  O 

•  • 

o  o 

'I  1  1  1  1  I 


Figure  20.  Measured  (Brum  et  al.  )  and  Predicted  Centerline 
Mean  and  rms  Axial  Velocity  for  UCI  CBCC . 
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Figure  21.  Measured  (Lightman  et  al.  )  and  Pri_1icted  Centerline 
Mean  and  rms  Axial  Velocity  for  APL  CBCC .  Air  Flow 
2  kg/s  and  Zero  CO-  Flow. 
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Figure  22.  Measured  (Lightman  et  al.  )  and  Predated  Centerline 
Mean  and  rms  Axial  Velocity  for  A.PL  CBCC .  Air  Flow 
2  kg/s  and  C0_  Flow  6  kg/hr. 


Figure  23.  Measured  (Lightman  et  al.2)  and  Predicted  Centerlin 
Mean  and  rms  Axial  Velocity  for  APL  CBCC.  Air  Flow 
2  kg/s  and  C0_  Flow  16  kg/hr. 


Figure  22  shows  the  FSP  to  occur  at  z/D  =  0.28  which  more 
closely  agrees  with  SS  and  KP  in  Figure  15  than  did  the  earlier 
experimental  data.  A  crucial  difference  in  the  two  measurements 
which  contributes  to  the  appreciable  difference  in  the  observed 
FSP  location  is  the  central-jet  exit  configuration.  Recent 
experiments  employed  a  well-designed  nozzle,  thereby  ensuring  a 
uniform  Wp.  Earlier  experiments  involved  a  straight  tube  15 
diameters  in  length  upstream  of  the  exit  plane,  resulting  in  a 
nonuniform  Wp.  while  this  difference  accounts  for  the  decrease 
in  the  FSP  location,  the  reason  for  the  increase  in  the  RSP 
location  is  not  clear. 

Pred ic t ions .  In  terms  of  the  present  predictions,  a  number 
of  observations  can  be  made.  Calculations  with  c M  and  e  modifi¬ 
cations  are  clearly  inferior  for  all  three  CO2  flows  in  all 
respects.  The  fact  that  the  three  predictions  do  not  differ 
significantly  for  16  kg/hr  reflects  that  the  central  jet  essen¬ 
tially  breaks  through  the  recirculation  region  and  does  not 
suffer  the  large  streamline  curvature  effects  associated  with  the 
recirculating  flow.  At  the  lower  CO2  flows,  the  dissipation 
modification  significantly  decreases  the  rms  velocity.  For  the 
mean  velocity,  the  zero  CO2  case  shows  good  agreement  for  the 
RSP  and  very  poor  agreement  for  th ?  velocity  profile  in  the 
reverse- t low  region.  For  the  6  kg/hr  CO2  case,  while  the  magni¬ 
tude  of  Wm  agrees  reasonably  well  with  the  measured  value,  its 
location,  as  well  as  the  locations  of  FSP  and  RSP,  are  all  con¬ 
siderably  overpred ic ted . 

The  differences  in  the  predictions  based  upon  standard 
cu  and  curvature-corrected  cu  are  not  clear  cut.  At  16  kg/hr, 
the  uncorrected  prediction  for  W  is  superior,  while  the  correction 
yields  better  w  prediction,  especially  near  the  RSP  (as  may  be 
anticipated).  At  6  kg/hr,  the  correction  again  predicts  w 
better  near  the  RSP.  For  W,  except  for  the  location  of  the  RSP, 
the  curvature  correction  appears  to  show  better  agreement  v  i th 
the  measurement  than  the  standard  model.  Both  calculations 
underpredict  the  location  of  Wm ;  the  magnitude  is  underpredicted 
by  the  standard  model  and  overpredicted  by  the  curvature 
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correction.  For  the  zero  CO2  flow  case,  the  uncorrected  predic¬ 
tion  appears  to  be  better  overall. 

CO2  Concentrations.  Comparison  of  the  predictions  with  the 
measurements  (Bradley  et  al.)30  of  CO2  mole  fractions  seen  in 
Figures  24  and  25  leads  to  similar  conclusions  regarding  the 
calculations.  The  prediction  of  the  standard  model  is  clearly 
superior  at  the  iiigher  flow  rate.  4t  6  kg/hr,  the  standard  model 
and  the  curvature  correction  do  not  appear  to  differ  appreciably 
from  each  other.  The  predictions  with  the  modifications  of 
c  u  and  r,  show  much  slower  decay  of  CO2  centerline  concentration 
than  do  the  measurements  and  other  predictions.  This  behavior  is 
consistent  with  the  significant  underprediction  of  w  and  the 
slower  decay  of  W  (noted  especially  for  6  kg/hr  in  Figure  22). 

d .  Cone  1 usions 

The  numerical  computations  based  upon  the  Reynolds-averaged 
equations  and  k-e  turbulence  model  have  demonstrated  that  the 
CBCC  represents  a  highly  complex  turbulent  flow.  The  predictions 
with  the  standard  k- s  model,  as  well  as  ad  hoc  modifications  to 
account  for  the  streamline  curvature  and  the  preferential  influ¬ 
ence  of  normal  stresses  on  dissipation  are  compared  with  the 
measurements  in  a  large-scale  and  a  small-scale  CBCC.  Signifi¬ 
cant  differences  between  the  predictions  are  noted,  especially 
for  the  dissipation-equation  modification  which  results  in  much 
greater  dissipation  in  the  large-scale  CBCC.  The  introduction  of 
the  curvature-corrected  cu  results  in  better  predictions  for 
certain  flowfield  regions  but  exhibits  a  dependence  on  the 
annular  and  central  flow  rates. 

3.  VO  R  T  BX  -  C  K  NT  hi  R  CH  ARACTK  RI  ST  I CS 

The  time-averaged  characteristics  of  the  vortex  center  in 
the  n*  ir-wake  recirculation  region  in  the  CBCC  by  finite- 
difference  computations  are  examined.  The  numerical  predictions 
under  isothermal  conditions  examine  the  influence  of  the 
turbulence-model  corrections  under  varying  annul  r  air  and 
central  CO 2  flows  on  the  axial  and  radial  locations  of  the  vortex 
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Figure  32.  Axial  Coordinate  of  Vortex  Center  end  Minimum  Cen¬ 
terline  Velocity  Variation  with  CO2  Flow  (APL  CBCC ) 
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TABLE  4 

FLOWFIELD  CONDITIONS 


UCI 

APL 

.  2 

-4 

2 

-4 

mAxlO 

A 

WA 

Re. xlO 

A 

ITKXlO 

A 

WA 

Re .  x  1 0 

A 

1.24 

7.78 

1.07 

100 

23.38 

17.8 

1.86 

11.67 

1.60 

200 

46.76 

35.7 

2.48 

15.57 

2.13 

300 

70.14 

53.5 

3.10 

19.46 

2.67 

400 

93.52 

71.3 

3.72 

23.35 

3.20 

500 

116.9 

89.1 

r— 

J 

3 

.  5 

-3 

(mF/mA )xlO 

VWA 

ReF/ReA 

m,,xlO 

r 

WF 

Re^xlO 

0.00 

0.00 

0.00 

0 

0.0 

0.0 

0.80 

0.53 

0.05 

in 

34.1 

16.3 

1.59 

1.06 

0.10 

167 

51.2 

24.4 

2.  39 

1.59 

0.15 

222 

68.2 

32.5 

3.20 

2.13 

0.20 

278 

85.3 

40.6 

3.98 

2.66 

0.25 

333 

102.4 

48.8 

4.78 

3.19 

0.30 

389 

119.4 

56.9 

5.58 

3.72 

0.35 

444 

136.5 

65.0 

6.37 

4.25 

0.40 

500 

153.6 

73.2 

556 

170.6 

81.3 

611 

187.7 

89.4 

streamline-curvature  corrections,  and  dissipation  modifications 
[see  Paragraph  II.2.b(2)].  The  range  of  the  annular  and  central 
flow  rates  investigated  for  both  the  large-  and  small-scale 
CBCC's  corresponds  to  all  the  flow  regimes  depicted  in  Figure  26. 
Table  4  lists  the  various  characteristics  of  the  flowfield  for 
both  the  APL  and  UCI  CBCC.  The  Reynolds  number  for  the  air  flow 
Re^  is  based  on  D  and  the  reference  velocity  Wd  in  the  duct;  the 
d  and  Wp  (See  Figure  13).  Variation  in  the  C02  flow  conditions 
is  due  to  the  input  parameters  which  were  varied  in  the  actual 
experiments  (i.e.,  for  UCI,  the  ratio  rap/m^  was  varied,  whereas 
for  APL,  flip  was  varied). 

b .  Results _ a n d  D i sc u s s i on 

(1)  APL  Configuration 

Figure  31  shows  the  relation  between  the  normalized  radial 
coordinate  of  the  vortex  center  and  CO2  flow  rate.  In  what 
follows,  the  turbulence  model  with  curvature  correction  is 
denoted  as  PC  and  the  one  with  both  curvature  correction  and 
dissipation  modification  is  denoted  as  PCE.  For  both  the  PC  and 
PCE  models,  the  radial  position  stays  approximately  fixed  at  a 
value  around  0.36.  This  is  consistent  with  the  preliminary 
results  in  Paragraph  l[.3.a(l)  (which  noted  that  for  CO 2  flow 
rates  of  0  to  8  kg/hr  and  air  flow  rate  of  2  kg/s  the  normalized 
radial  coordinate  is  a  constant  of  approximately  0.35).  At  a 
flow  rate  near  11  kg/hr,  however,  both  models  exhibit  a  slight 
increase  in  distance  from  the  centerline.  This  move  could  be  a 
result  of  the  growing  prominence  of  the  central  jet  at  CO2  rates 
of  12  to  16  kg/hr.  (Reverse  mass  flow  on  the  centerline 
disappears  at  16  kg/hr  with  an  annular  flow  rate  of  2  kg/s.) 

The  slight  movement  of  the  radial  coordinate  is 
insignificant  in  comparison  to  the  migration  of  the  axial 
coordinate  zc,  shown  in  Figure  32.  In  the  region  of  cente-line 
reverse  flow,  the  axial  coordinate  moves  gradually  downstream 
with  increasing  CO 2  flow,  then  migrates  back  d  the 

centerbody.  The  peak  value  of  zc  occurs  for  .1  ;  iow  rate  ot 
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of  zero  central  flow  is  0.3  in  both  calculations.  At  the  highest 
CO2  flow  rate,  the  FREP  result  is  about  20  percent  higher  than 
the  TEACH  result.  Since  the  TEACH  calculations  are  based  upon 
the  physically  more  realistic  k- e  turbulence  model  than  the 
constant  eddy  viscosity  model  used  in  the  FREP  calculations, 
present  results  for  the  axial  location  of  the  vortex  center  may 
be  more  accurate.  The  identical  value  for  the  radial  location, 
however,  implies  that  the  details  of  turbulence  model  do  not  have 
much  effect  on  how  far  the  vortex  center  is  displaced  off  the 
centerline.  This  conclusion  is  consistent  with  the  observation 
of  Ko  and  Chan^G  that  the  radial  position  of  the  vortex  center  in 
their  unconfined  annular  jet  study  was  essentially  independent  of 
the  momentum  flux  of  the  annular  jet  and  hence  the  pressure 
available  for  the  entrainment  behind  the  centerbody  face.  For 
further  implications  of  our  results  vis-a-vis  those  of  Ko  and 
Chan, 36  we  recall  our  discussion  in  Reference  5. 

( 2 )  Objectives  of  Present  Study 

The  foregoing  results  were  limited  to  the  APL  CBCC  and  only 
a  small  range  of  the  central  CO2  flow  rates  for  a  fixed  annular 
air  flow  rate  was  considered.  Although  these  calculations  dealt 
with  isothermal  flowfields,  a  detailed  examination  of  the  vortex 
center  characteristics  could  offer  additional  insights  into  the 
phenomenon  of  flame  stabilization.  Accordingly,  the  main  objec¬ 
tive  of  the  present  study  is  a  parametric  investigation  of  the 
vortex  center  for  varying  air  and  CO2  flow  rates  under  different 
turbulence-model  corrections  in  both  APL  and  UCI  CBCC's.  A 
second  objective  is  to  examine  the  influence  of  the  variations  in 
flow  rates,  turbulence-model  parameters,  and  combustor  scaling  on 
the  magnitude  and  location  of  the  minimum  (which  is  negative  when 
centerline  flow  reversal  is  present)  centerline  axial  velocity. 

(  3  )  S_cope 

This  study  deals  with  the  time-averaged  flowfield  pre¬ 
dictions  and  is  based  on  certain  modifications  to  the  standard 
TEACH  numerics.  These  include  a  power-law  differencing  scheme, 
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flux  represents  the  positive  axial  flux  contained  between  the 
zero  velocity  contour  and  the  time-averaged  separated  streamline. 
In  other  words,  the  radially  outward  flux  cannot  cross  the 
time-averaged  separated  streamline  into  the  annular  main  flow. 

The  distribution  of  the  turbulence  kinetic  energy  on  the 
zero  velocity  contour  is  shown  in  Figure  28.  It  is  seen  that  the 
local  maximum  of  the  kinetic  energy  on  the  contour  within  the 
recirculation  region  occurs  at  the  vicinity  of  the  vortex  center. 
An  examination  of  the  radial  distribution  of  the  kinetic  energy 
at  five  axial  locations  shown  in  Figure  29  indicates  that  except 
for  the  axial  location  closest  to  the  centerbody  face,  the  tur¬ 
bulence  kinetic  energy  peaks  near  the  zero  velocity  contour. 

These  observations  appear  to  indicate  that  the  time-averaged  tur¬ 
bulence  activity  is  highest  along  the  zero  axial  velocity  contour 
inside  the  recirculation  region.  It  is  physically  reasonable  to 
expect  that  the  region  of  highest  turbulence  activity  would 
represent  the  region  of  optimum  mixing.  Although  the  numerical 
computations  reported  here  correspond  to  isothermal  flowfields 
and  do  not  necessarily  apply  to  reacting  flowfields,  our  current 
viewpoint  has  interesting  implications  for  bluff-body  flame 
s tabi 1 i zat ion . 3 3- 35  it  is  tempting  to  associate  the  vortex 
center  with  a  local  hot  spot  (and  temperature  peak)  and  the  zero 
velocity  contour  with  the  optimal  region  of  flame  propagation. 

T ime- Averaged  Vortex  Center.  The  radial  and  axial  coor¬ 
dinates  (normalized  with  respect  to  the  centerbody  diameter)  of 
the  vortex  center  for  different  central  CO2  flow  rates  are  shown 
in  Figure  30.  We  note  that  the  radial  location  of  the  vortex 
center  is  invariant  with  CO2  flow  rate.  The  axial  coordinate, 
however,  is  seen  to  increase  slowly  with  CO2  flow  rate.  A  com¬ 
parison  of  these  results  with  those  of  our  previous  study1’  with 
the  FREP  code  (which  employed  a  constant  value  of  turbulent  eddy 
viscosity)  shows  the  somewhat  intriguing  result  that  the  pre¬ 
dicted  locations  of  the  vortex  center  in  both  TEACH  and  FREP 
modeling  are  not  much  different.  Indeed,  the  r^lial  location  is 
predicted  as  0.35  by  both  codes.  The  axial  location  for  the  case 
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center,  as  well  as  on  the  magnitude  and  location  of  the  minimum 
centerline  axial  velocity.  These  results  emphasize  and  clarify 
the  complex  flowfield  interactions  of  the  annular  and  central 
streams  in  the  near-wake  region. 

a .  Recirculation-Zone  Structure 

The  nature  of  the  toroidal  recirculating  vortex  in  the  CBCC 
is  illustrated  in  Figure  26  which  shows  the  computed  velocity 
field  in  the  APL  CBCC.  The  computations  involve  annular  air 
flow  and  central  CO2  flow,  with  the  figure  depicting  the  flow- 
fields  corresponding  to  (a)  no  central  flow,  (b)  small  central 
flow,  and  (c)  large  central  flow.  The  recirculation  region 
comprises  two  toroidal  vortices:  a  primary  vortex  and  a  smaller 
secondary  vortex  near  the  central  jet.  The  relative  sizes  of 
these  vortices  and  the  locations  of  the  vortex  centers  are 
dependent  upon  the  annular  and  central  flow  rates  and  the 
resulting  degree  of  mixing  and  entrainment  present  in  the  system. 
At  zero  or  small  central  flow,  the  annular  flow  is  dominant, 
entraining  the  CO2  flow  and  causing  a  flow  reversal  along  the 
centerline.  As  the  strength  of  the  central  jet  increases,  it 
gradually  overcomes  the  centerline  reverse  flow,  and  the 
secondary  vortex  grows  relative  to  the  primary  vortex.  Eventu¬ 
ally  the  central  jet  becomes  strong  enough  to  eliminate  the 
centerline  reverse  flow  entirely  and  starts  to  entrain  the  an- 
ular  flow.  Off  centerline  reverse  flow  and  primary  vortex  are 
still  present,  the  latter  at  diminished  strength  however. 

( 1 )  Preliminary  Studies 

Associated  with  the  primary  vortex  is  the  time-averaged 
contour  of  zero  mean  axial  velocity.  Figure  27  presents  this 
contour  and  the  mean  radial  velocity  (U)  distribution  thereon  for 
the  case  of  zero  central  flow.  The  point  of  zero  radial  velocity 
on  the  contour  is  the  vortex  center.  As  may  be  expected  or  phy¬ 
sical  grounds,  the  flow  is  radially  outward  to  the  left  and 
radially  inward  to  the  right  of  the  vortex  center  on  the  zero 
velocity  contour.  Note,  however,  that  the  radially  outward  mass 
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Figure  25.  Measured  (Bradley  et  al.^O)  and  Predicted  Centerline 
CC>2  Mole  Fraction  for  APL  CBCC .  Air  Flow  2  kg/s  and 
CO2  Flow  15  kg/hr. 
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Measured  (Bradley  et  al.30)  and  Predicted  Centerline 
CO2  Mole  Fraction  for  APL  CBCC.  Air  Flow  2  kg/s  and 
CO2  Flow  f<  kg/hr. 


7  kg/hr  for  PCE  and  10  kg/hr  for  PC.  This  result  can  be 
explained  by  the  competition  between  the  annular  and  central 
flows.  As  the  dominance  of  the  central  flow  increases,  the 
entrainment  by  the  annular  flow  gradually  diminishes  and 
eventually  the  annular  flow  begins  to  get  entrained  by  the 
central  flow  rather  than  the  converse,  which  is  true  at  low  CO2 
flows.  It  is  here  that  zc  peaks  and  the  primary  vortex  (and  its 
center)  is  pushed  downstream.  When  entrainment  of  the  annular 
flow  begins,  the  primary  vortex  is  pulled  toward  the  centerbody. 
Near  the  point  at  which  negative  centerline  flow  disappears 
(around  16  kg/hr),  a  secondary  peak  is  also  found  for  zc.  it 
appears  that  the  detachment  of  the  primary  vortex  from  the 
centerline  promotes  a  small  downstream  motion  of  the  vortex 
center.  After  this  peak  is  reached,  all  increases  in  CO2  flow 
draw  the  primary  vortex  center  toward  the  centerbody. 

While  the  general  trend  of  the  path  of  the  vortex  center  is 
consistent  between  the  PC  and  PCE  models,  the  actual  magnitude  of 
the  axial  coordinates  and  the  precise  location  of  the  peaks  vary. 
The  PCE  model  tends  to  predict  values  further  downstream  than  the 
PC  model.  This  same  trend  is  evident  in  the  location  of  tne 
minimum  centerline  velocity  (zm).  This  location  steadily  moves 
downstream  with  increasing  CO2  flow,  although  the  actual  magni¬ 
tude  of  zm  is  dependent  upon  the  presence  of  the  dissipation 
modification.  According  to  Paragraph  II.2.c(3),  the  PC  model 
tends  to  predict  the  minimum  velocity  location  and  magnitude  more 
accurately  for  the  APL  CBCC. 

If  air  flow  rate  is  varied,  holding  CO2  flow  rate  constant, 
the  radial  coordinate  of  the  vortex  center  again  remains  fixed 
near  0.36  at  CO2  flow  rates  of  0  and  6  kg/hr  (see  Figure  33). 

This  result  is  consistent  with  the  findings  of  Ko  and  Chan36  that 
in  their  unconfined  annular  jet  the  radial  position  of  the  vortex 
center  is  not  dependent  on  the  annular  momentum  flux.  In  fact, 
the  radial  coordinate  normalized  by  the  duct  diameter  is  0.19  and 
this  compares  very  well  with  the  value  of  0.18  reported  by  Ko  and 
Chan.  Their  results  for  no  central  flow  and  increasing  annular 
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flow,  however,  show  a  shift  of  the  vortex  center  toward  the 
centerbody.  This  study,  on  the  other  hand,  demonstrates  no 
change  in  the  axial  position  of  the  vortex  center  with  increasing 
air  flow  and  no  CO2  flow.  The  differences  may  be  explained  by 
noting  that  Ko  and  Chan  used  an  unconfined  jet.  Such  a 
configuration  allows  for  ambient  air  entrainment  at  the  annular 
jet,  something  which  a  ducted  CBCC  cannot  provide.  With  the 
addition  of  the  central  jet,  however,  the  extra  potential  for 
entrainment  is  present,  thus  giving  results  similar  to  those  of 
Ko  and  Chan  (i.e.,  a  movement  of  the  axial  coordinate  toward  the 
centerbody  with  increased  annular  flow).  At  lower  annular  flows, 
entrainment  of  the  annular  jet  by  the  central  jet  is  a  more 
prominent  feature,  thus  causing  the  vortex  center  to  be  located 
upstream  of  the  vortex  center  at  higher  air  flows. 

( 2 )  UCI  Conf igurat ion 

Figures  34  and  36  reinforce  the  APL  results  [discussed  in 
II. 3. b( 1)]  in  that  the  radial  coordinate  of  the  vortex  center  is 
essentially  fixed  at  approximately  0.36,  regardless  of  CO2  and 
air  flow  rates.  Once  again,  a  slight  increase  in  distance  from 
the  centerline  is  observed  as  the  CO2  flow  rate  approaches  the 
point  of  elimination  of  negative  centerline  velocity.  This  shift 
could  be  caused  by  the  elongation  and  lifting  of  the  primary  vor¬ 
tex  from  the  centerline. 

Figure  35  shows  that  the  variation  of  the  axial  coordinates 
as  predicted  by  the  PCE  model  follows  the  same  trend  as  in  the 
APL  case.  The  first  peak  occurs  consistently  (for  both  APL  and 
UCI  PCE  models)  at  a  Wp/W^  ratio  of  approximately  1.3.  The 
secondary  peak  occurs  at  a  Wp/W^  ratio  of  approximately  2.6. 

The  uniformity  at  this  point  is  not  as  strong,  however.  The 
general  trend  for  the  variance  of  zc  with  ftp  is  that  with 
increasing  air  flow,  the  peak  values  of  zc  are  less  pronounced, 
due  to  the  increased  prominence  of  the  annular  jet.  Results 
from  the  PC  model  (not  shown)  concur  with  those  of  the  PCE 
model.  The  values  for  the  location  of  the  mini.],  urn  centerline 
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Figure  36.  Radial  Coordinate  of  Vortex  Center  Variation  with  A 
Flow  (UCI  CRCC )  -  PC  Model. 


velocity,  also  shown  in  Figure  34,  confirm  the  APL  results  in 
that  zm  is  nearly  linearly  dependent  upon  the  CO2  flow  rate.  An 
increase  in  annular  flow,  as  expected,  pushes  the  zm  value  closer 
to  the  centerbody.  A  graph  of  zm/D  versus  Wf/Wa  (not  shown) 
would  reveal  that  the  results  for  all  the  three  air  flows  (7.8, 
15.6,  and  23.4  m/s)  essentially  coincide. 

Figure  37  presents  the  variation  of  the  axial  coordinate  of 
the  vortex  center  with  annular  air  flow.  The  zc  value  is 
essentially  constant  for  low  CO2  rates,  analogous  to  the  0  kg/hr 
CO2  flow  in  APL  (see  Figure  32),  but  at  a  higher  value  of  0.38 
compared  to  0.36  for  APL.  At  higher  CO2  flow  rates,  the  trends 
are  not  similar  to  those  of  the  APL  CBCC.  One  possible  explana¬ 
tion  is  the  difference  in  the  Wp/WA  ratios.  In  the  APL  con¬ 
figuration,  the  ratios  ranged  from  0.44  to  0.73  in  the  region 
which  exhibited  a  decrease  in  zc  with  increased  air  flow.  The 
one  point  which  appeared  to  be  dominated  by  central  flow  had  a 
Wp/WA  of  2.19.  All  of  the  points  of  the  UCI  curves  at  the  three 
higher  CO2  flows  have  the  Wp/WA  values  of  2.13  or  higher. 

( 3 )  Comparative  Study 

Figures  38  and  39  demonstrate  the  correlation  between  the 
APL  and  three  CO2  flow  rates  of  UCI.  The  normalized  minimum  cen¬ 
terline  velocity  as  a  function  of  normalized  central  flow  has, 
for  the  most  part,  a  constant  relationship  in  all  runs  (see 
Figure  38).  The  ratio  between  zm  and  zc,  shown  in  Figure  39, 
appears  consistent  among  all  runs  when  compared  to  the  ratio  of 
the  inlet  flow  velocities,  although  these  results  do  not  coincide 
as  well  as  those  in  Figure  38.  Similar  trends  can  be  seen  be¬ 
tween  the  PC  and  PCE  models  in  both  correlations,  though  for  each 
case  only  one  model  is  shown.  This  relation  between  the  minimum 
centerline  velocity  and  vortex  center  location  could  prove  help¬ 
ful  in  experimental  techniques  where  inlet  flow  rates  and  cen¬ 
terline  velocities  are  easier  to  measure  than  is  the  vortex 
center  location. 

The  results  from  this  study  extend  the  preliminary  results 
in  Figure  30  and,  as  noted  in  Paragraph  Tl.3.a.(l),  appear  to 
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Figure  37.  Axial  Coordinate  of  Vortex  Center  Variation  with  Air 
Flow  (IJCI  CBCC )  -  PC  Model. 
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Relationship  of  Axial  Location  of  Minimum  Centerl 
Velocity  and  Central  Jet  Flow  (APL  and  UCI  CBCC ) 
PCS  Model. 


coincide  with  those  of  a  previous  study  on  primary  vortex  charac¬ 
teristics  by  Krishna mu r thy. ^  Using  the  Field  Relaxation  Elliptic 
Procedure  (FREP)  Code,  Kr ishnamurthy  obtained  the  constant  value 
of  0.35  for  the  normalized  radial  coordinate  at  various  central 
jet  flows.  Also,  the  axial  coordinate  of  the  vortex  center 
moved  downstream  with  increased  central  flow  in  the  range 
observed  (0  to  8  kg/hr  CO2  flow  and  2  kg/s  air  flow).  In  a  later 
investigation,  Kr ishnamurthy ,  et  al.3  compared  the  results  of  the 
standard  TEACH  code  to  the  experimental  results  of  Lightman,  et 
al.2  for  annular  flow  of  2  kg/s  and  no  central  jet.  The  radial 
again  was  found  to  be  consistent  at  a  r/D  value  of  0.35,  in 
agreement  with  present  results.  For  the  axial  coordinate, 
however,  the  standard  TEACH  code  considerably  underpredicted  the 
location  of  the  vortex  center  as  compared  to  the  experimental 
data  (see  Table  5). 

TABLE  5.  AXIAL  COORDINATE  OF  VORTEX  CENTER  IN  VARIOUS  MODELS 

( APL  CBCC) . 


Mode  1 

Normalized  Axial  Location 

Experiment2 

0.43 

FREP5 

0.30 

TEACH  (standard)3 

0.26 

TEACH-PC 

0.31 

TEACH-PCE 

0.31 

All  of  the  previous  and  present  studies  have  underpredicted 
the  axial  location  of  the  vortex  center  to  some  degree. 
Streamline-curvature  correction  and  dissipation-equation  modifi¬ 
cation  have  provided  the  closest  approximations  thus  far  to  the 
experimental  results.  With  the  central  ;et  added,  the  general 
trend  has  been  that  the  PCE  Model  predicts  a  vortex  center  loca¬ 
tion  further  downstream  than  the  PC  model.  This  indicates  that 
better  agreement  with  experimental  results  for  increased  CO2 
flow  is  likely  with  the  PCE  model. 


Conclus ions 


c . 

The  radial  coordinate  of  the  primary  vortex  center  remains 
essentially  constant  regardless  of  air  or  CO2  flow  rates.  The 
axial  coordinate  moves  downstream  with  increasing  central  flow 
and  constant  annular  flow  until  entrainment  of  air  by  CO2  begins 
then  zc  migrates  toward  the  centerbody.  With  changing  air  flow 
and  zero  CO2  flow,  the  axial  coordinate  stays  constant.  The 
addition  of  central  flow  causes  zc  to  shift  toward  the  center- 
body,  if  the  annular  jet  is  dominant.  The  differences  between 
the  UCI  and  APL  combustors  are  seen  mostly  in  the  magnitude  of 
the  variables;  normalization  shows  the  same  overall  trends 
preserved  between  the  two,  although  not  exactly. 
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SECTION  III 


TIME-DEPENDENT  COMPUTATIONS 

This  section  deals  with  the  calculation  of  the  CBCC  flowfield 
by  the  solution  of  the  time-dependent  Mavier-Stokes  equations 
using  the  AFWAL/FIMM  unsteady  code. 


1.  INTRODUCTION 

Proper  specification  of  boundary  conditions  is  critical  to 
the  numerical  solution  of  subsonic  flowfields.  This  section 
presents  a  systematic  study  of  the  influence  of  various  inflow 
and  outflow  boundary  conditions  for  the  numerical  solution  of  the 
time-dependent  Navier-Stokes  equations.  Combinations  of  several 
representative  inflow  and  outflow  boundary  conditions  are  applied 
to  obtain  the  solution  of  subsonic  unsteady  flow  in  a  centerbody 
combustion  chamber  of  finite  computational  domain.  The  absence 
of  solutions  of  nonlinear  boundary  value  problems  of  this  class 
is  explained  both  by  the  difficulties  in  constructing  models  of 
unsteady  subsonic  flows  of  a  viscous  gas  in  a  finite  region  (due 
to  the  specific  requirements  imposed  on  the  boundary  conditions), 
and  by  the  problem  of  correct  resolution  of  the  acoustic 
processes  by  the  difference  scheme.  For  subsonic  flow,  the 
propagation  of  disturbances  upstream  against  the  flow  direction 
complicates  the  situation  significantly  and  requires  more  care  in 
the  choice  of  inflow  and  outflow  boundary  conditions.  Following 
the  approach  of  Reference  37  the  acoustic  characteristics  of 
various  boundary  conditions  are  investigated,  which  enables  one 
to  study  their  effect  on  the  longitudinal  oscillations  in  the 
domain  of  interest.  In  contrast  to  subsonic  flow,  supersonic 
flow  leaves  little  choice  in  the  selection  of  the  boundary 
conditions  and  all  physical  variables  should  be  specified  at  the 
inflow  boundary  and  none  at  the  outflow  boundary. 

At  a  subsonic  inflow  hounu  ry,  the*  Pav ier-  '  tokes  equations 
require  the  specification  of  three  boundary  conditions  and  the 


subsonic  outflow  boundary  requires  the  specification  of  only  one 
boundary  condition.  Therefore,  in  the  numerical  solution  of  the 
problem  the  variables  at  the  boundary  points  must  be  calculated 
by  introducing  additional  difference  algorithms,  containing  to 
some  degree  elements  of  spatial  extrapolation  from  the  interior 
points  of  the  domain.  Various  ex t rapola t ion  procedures  of  the 
type  (  3 <{>  /9  n  )  n  +  ^  =  0  or  (  3  ^  /a  n  2)  n  +  1  _  g  can  be  use(j  to  compute 

the  additional  unspecified  physical  variables  required  by  the 
finite-difference  scheme  at  the  boundaries  (here  n  is  in  the 
direction  normal  to  the  the  boundary  and  n  corresponds  to  the  the 
time  step).  However,  there  is  no  rigorous  mathematical  theory 
available  to  justify  such  extrapolation  procedures. 

2.  BOUNDARY  CONDITIONS 

In  what  follows  various  inflow  boundary  conditions  and  their 
finite-difference  analogues  are  described  and  the  boundary  con¬ 
ditions  are  applied  at  both  the  predictor  and  corrector  steps  of 
MocCormack's  finite-difference  algorithm.  The  conditions  are 
described  only  for  the  (n+1)  corrector  step.  The  procedure  for 
the  (n+1/2)  predictor  step  is  identical  to  the  (n+1)  corrector 
step,  with  the  time  step  superscripts  (n+1;  n+1/2)  replaced  by 
(n+1/2;  n).  The  subscripts  denote  the  x  (axial)  and  r  (radial) 
grid  indices  and  KL  is  the  last  grid  index  in  the  x  direction. 

a  .  Inf  low  Jjounda  ry  Conditions 

The  following  paragraphs  discuss  the  various  combinations  of 
the  boundary  conditions  at  the  inlet. 

( 1 )  Boundary  Conditions  due _ to  Serra 

Serra^  suggests  the  use  of  total  pressure,  total 
temperature,  and  fLow  angle  at  the  inflow  boundary  for  successful 
modeling  of  internal  gas  flows. 

l’t  ,  Tfc  ,  V  =  0  and  =  0  (8) 


79 


AD-A153  682 
UNCLASSIFIED 


TINE-AVERAGED  AND  TINE-DEPENDENT  CONFUTATIONS  OF 
ISOTHERNAL  FLOUF I  ELDS  IN.  .  <U>  DAVTON  UNIV  OH  RESEARCH 
INST  L  KRISHNAHURTHV  ET  AL.  DEC  84  UDR-TR-84-88 
AFHAL-TR-84-2881  F23615-82-K-2252  F/G  28/4 


The  f i n i te-d i f ference  formulation  of  Eq .  (8)  is  given  by 


Tn  +■  1 

1,1 

Tt  - 

(Y-l)  [u!?  +  Jl  2/(  2y  R ) 

(9) 

pn+l  _ 
l,j 

Pt/ 

f  1+  (y-1)  [u"+Jl  2/(2YRT"fh}  y/y_1 

1  r  J  1  r  J 

(  10  ) 

v-1  = 
1  ,  1 

0.0 

(  11  ) 

ll"  +  1  = 

1  ,  3 

u"  +  1 
2,1 

(  12  ) 

n  +  1  - 

PM  - 

pn  +  l 
1,  j 

/( RTa+J  ) 

f  J 

(13) 

Knowing 

u?+1. 
1 ,  1 

and  from  Eqs.  (11)  and  (12)  is 

1,1  1,1 

determined 

from 

Eq.  (9)  and  P^+  4  ,  +  4  follow  from  Eqs. 

1,1  1,1 

(10) 

and  (13),  respectively. 


(2)  A  Minor  Variant  of  Serra's  Conditions 

A  slight  variant  of  Serra's  boundary  conditions  results  if 
the  condition  3  U/3  x  =  0  is  replaced  by  3pU/3x  =  0.  We  have 

Pt'  Tt'  v  =  °'  and  fir  =  0  (14) 

The  f i n i te-d i f f e rence  formulation  of  Eq.  (14)  is  the  same  as 
Eg.  (9)  to  (13)  except  for  Eq .  (12)  which  is  replaced  by 


(  a  )  _Ei_g  o  r  ous  Imp  lemon  t  at  tin 

Because  of  the  nonlinearity  in  Eq.  (15),  Kqs.  ( 9  )  ,  (  1 0  )  ,  (  1  3  ) 
and  (15)  form  a  set  of  coupled  nonlinear  equations  whose  solution 
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procedure  is  facilitated  by  the  following  method: 


oJtj  -  KlTUj/R 

(  16) 

vT,i  -  K/<] 

(17) 

K2ITu5'2/T  ^  TUj  *  Tt 

(  18) 

where 

and 

K2  =  [  ( y -1 )  R/2y  1  [  K/K :  1  2  • 

The  solution  of  Eq.  (18)  is  accomplished  iteratively  through 
quas i 1 inearization  as  seen  below: 


(Tt-  6K2/[Tj;>q  )  /(  1-5K2/[T"^]^  ) 


(19) 
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By  making  use  o£  the  predictor  step  values  of  U  and  V  Eqs. 

(20)  and  (21)  can  be  solved  for  and  ,  independent  of 

1 1  1  w  J 

the  updated  values  u’14^  and  V1?4^  obtained  from  Eqs.  (22) 

1  /  1  f  »  3 

and  (23).  Density  is  obtained  from  Eq .  (24).  Though  this 

approach  reduces  the  program  complexity,  it  is  found  to  produce 
an  entirely  different  solution  from  that  due  to  the  rigorous 
implementation  of  the  inlet  conditions  as  outlined  by  Eqs.  (15  to 
19)  [Bee  paragraph  T[[.2.a(2a)]. 

(3)  Boundary  Conditions  Due  to  Oliger  and  Sundstrom 

Oliger  and  Sundstrom^0-41  have  shown  that  the  specification 
of  p,  H,  and  V  at  the  inflow  is  well-posed  for  the  linearized 
system  of  equations.  We  have 

p ,  U,  V  =0  and  |I  =  0.  (25) 

O  X 

The  finite-difference;  analog  of  Eq.  (25)  is  given  by 


n+  1 

PU] 

-  f  (  t) 

(  26) 

H?41 

1  .  J 

=  g  ( t ) 

(  27) 

i  f  i 

II 

O 

• 

o 

(  28) 

1  1  »  j 

2f  J 

(  29) 

P?  +  1 

1  f  1 

n+1  nfl 

-  ^  1  ,  j  T  L ,  j 

(  30) 

where  f(t)  and  g(t)  are  given  functions,  which  are  fairly  smooth 
and  monotonic  for  t  e  ( t  ^ ,  t2>  and  are  practically  unchanged 
within  a  characteristic  period  of  oscillation  t}  =  2tt /w  <  t2~t^. 

( 4 )  Bound a r /  Conditions  Due  to  Fedorchenko 
37,42 


A.  T.  Kednrehenko 


successfully  modeled  unsteady  subsonic 


flows  of  a  viscous  case  with  the  inflow  boundary  conditions 
involving  specification  of  T,  pU,  and  pV: 


T,  p 

ii 

> 

a. 

°'  and  & 

■ 

(31) 

The 

f i n i te-d 

i f  f erence 

formulation  is  given  by 

Tn  +  i 

=  Cons 

tant 

(  32) 

1  /  1 

n  +  1 
P  i 

va  +  1  = 

0.0 

(  33) 

1- J 

1  /  1 

n+1 
P  1  r  j 

u"  +  1  = 
1.1 

f  (  t) 

(  34) 

Pa  +  1 

=  p"+1 

(35) 

1  /  1 

2.  J 

.  n+1 

Pl,j 

-  O 

(  36) 

b. 

Outflow 

Boundary 

Conditions 

In  what  follows  we  describe  different  variants  of  essentially 
what  amounts  to  the  specification  of  exit  pressure  (Pe)  as  the 
required  outflow  boundary  condition. 

( 1 )  Non-Reflecting  Condition 

Rudy  and  Strikwerda43  suggest  the  use  of  a  non-reflecting 
boundary  condition  which  is  based  on  the  outflow  characteristic 
of  the  Euler  equation  of  motion  for  application  to  viscous  gases. 
The  main  purpose  of  the  condition  is  to  prevent  upstream  propaga¬ 
tion  of  acoustic  waves  from  the  outflow  boundary  and  thereby 
accomplish  rapid  convergence  to  the  desired  steady-state  solution 
using  time-dependent  Navier-Stokes  equations.  Thus 


3  P  3  IJ 

Tt  pa  5T 


+  a(P-P  )  =  0, 


32U 


3  x 


32V 


3  x 


3  2T 


3  x 


=  0. 


(37) 


Because  of  the  time  dependence  of  the  non-reflecting  boundary 
condition  it  is  necessary  to  describe  the  finite-difference 
formulation  of  both  the  predictor  and  corrector  steps. 
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The  Predictor  Step: 


3  2  a  n  + 

'  3x2~ 


0 


(  38) 


where  <*>  =  U,  V,  and  T 

p  n  +- 1  /  2  _  pH  4-  p  n  (  y  RT 

KL,]  KL ,  j  pKL,j'Y 


KL,  j 


)l/2([Jn+l/2_  lJn 


KL,  j 


KL,  j 


•  ) 


(  39) 


-  a  A  T  (  P 


n 


The  Corrector  Step: 


a 2. n+1 

L*  i ta 

3x2 


0 


KL,  j 


Pe) 


(  40) 


where  4>  =  U,  V,  and  T 

+  i  .  =  0. 5(  Pnf'1//2+  pn  .) 

•  ^  v  -  r/  r  4  *•  r  r  r  4  ' 


KL,  j 


KL,  J 


KL ,  j 


+  °7l.?  1/2(,jkl!:-  uTi,i2 


(41  ) 


»4TI|,;I!f-pe) 


Instead  of  the  vanishing  second-order  derivatives  given  by 
Eqs.  (38)  and  (40),  Linear  extrapolation  of  <j>  =  U,  V,  and  T  is 
incorporated  at  the  exit,  as  seen  below: 

*KL,j  "  [(XKL,]  ~  XKL-l,  j)/(XKL~l,  j"  XKL- 2  ,  j  )  1  (  ♦  KL- 1  ,  j  “  <)>KL-2,j) 


KL- 1 , j 


(42) 


(2)  Constant  Pressure  Cond i t ion 

The  boundary  condition  most  widely  used  in  many  gasdynamic 
applications  involves  the  specification  of  constant  pressure. 


P  =  Pn, 


3  n 

iTx 


a  v 

;fx 


,n  t- 1 


3  T 
3~x 

,n  f  1 


=  0 


k i. ,  j 

.n  *-  1 
<1  ,  1 


e '  nKI.,]  UKI,-l,j' 


TnM 


KL-  1  ,  j 


a  n  1 1 


n  t- 1 
P  KL,  ] 


Vn  +  l.  =  v'14'1  ., 

KL,]  VKL- l  ,  ] 


~  n  4- 1 

kl,]  (<;;fJ> 
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(  44  ) 


(  3 )  Split  Boundary  Condition 


Hasen44  suggested  the  use  of  a  split  condition  for  pressure 
as  the  outflow  condition. 

( a )  Rigorous  Implementation 


P  =  PP, 


LE  =  n 

3  x  U ' 


r  <  R, 


R  <  r  <  R 
c  d 


3  U 
3  x 


3  V 

3  X 


3  T 
Tx 


=  0 


split  condition 


(45) 


where  Rc  and  R<j  are  the  radii  of  the  centerbody  and  duct  wall  of 
the  combustor,  respectively. 


P n  +  ^  =  P 

”kL , j 

n+1  =  pn  +  l 

KL,  j  KL- 1 , j 


if  r  <  R, 


if  R  <  r  <  R, 
c  d 


(46) 


n  +  1 


n  +  1 


A..  _  _  .  -  n  +  1  _  nn+l  ..--n+l  . 

^  KL ,  j  ^  KL- 1 , j  '  P  KL , j  PKL , j  RTKL , j 

where  <f>  =  U,  V,  and  T. 


( b )  Alternative  Implementation 

Instead  of  evaluating  density  from  the  equation  of  state,  Eq. 
(46),  Scott  and  Hankey^9  used  the  following  method: 


n  +  1 


-  =  pn+l  /  n+i 

Kt"j  ~  Kt-'j  (RTtCL-l,jll+Y(T'2,H  /2lKL-l,j} 


,n  +  l 


(47) 


Furthermore,  after  the  determination  of  the  primitive 
variables,  the  programming  of  Reference  39  in  the  evaluation  of 
the  conservative  variable,  pe,  is  as  follows: 


[pel 


n+1 
KL,  j 


=  p 


n  +  1 
K I 


,n  +  l 


n  +  1 


(  48  ) 
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The  implication  of  Eq .  (48)  is  that  more  energy  is  added  at 
the  outflow  boundary  by  an  amount  equal  to  l/?(u2  +  )  because 

specific  internal  energy  is  defined  as  e  =  CVT  +  l/>(tj2  +  )  . 

( 4 )  Relaxation  Boundary  Condition 

The  following  variant  involving  the  specification  of  pressure 
at  the  exit  is  used  in  Reference  45  for  the  modeling  of 
Kelvin-Helmholtz  instability. 


P 


n+1 
KL,  j 


P  +  {XKL“  XKL-l)(Pe“P)/Rd' 


where  P  is  obtained  from  the  linear  extrapolation;  and 


- — =  0,  where 
D  x  2 


U,  V,  and  T. 


(  49  ) 


(  5 )  Boundary  Conditions  For  Reverse-Flow  _Ex_i_t 

A.  T.  Fedorchenko42  suggests  the  use  of  specifying  V,  T,  and 
P  if  there  is  a  flow  reversal  at  the  exit. 


3  U  =  3  T  9  V 

e  '  3~x  3  x  3~x 


=  0 


a  n  d 


(  50) 


V  = 


0,  T  =  T  ,  P  =  P  and 
e  e 


3  p  U  _ 


n  +  1 


3  x 


=  o  if  u*1;  .  <  0 

KL,  j 


SoL_  id  -  VVa  1_  l_  _  Boundary  Cond  i  t  ions 


The  eenterbody  and  duct  walls  are  treated  as  no-slip, 
isothermal  surfaces,  involving  the  specification  of  Tw  and  Uw  = 

Vw  =  0.  The  remaining  thermodynamic  variables  Pw  and  pw  are  eval 
uated  from; 


P 

w 


Rp  T  , 
w  w 


3  P 
3  n 


0 


(51  ) 


or  ointinuif/  equation,  where  n  is  in  the  direction  normal  to  the 


d .  Centerline  Boundary  Conditions 

The  boundary  conditions  on  the  axis  of  symmetry  are  given  by 

IT  =  5T  =  3T  =  0  and  V  -  0  (52 

e .  Acoustic  Characteristics  of  Boundary  Conditions 

By  linearizing  the  inflow  and  outflow  boundary  conditions  it 
is  possible  to  estimate  the  acoustic  properties  of  the  longitu¬ 
dinal  oscillations  in  the  domain  of  interest.  Following  the 
aproach  of  A.  T.  Fedorchenko , 3 7  the  relative  value  of  normal  con¬ 
ductance  0  and  reflection  coefficient  3  can  be  estimated  for  the 
given  inlet  boundary  by  considering  small  adiabatic  fluctuations: 


U  =  U  +  U'  ,  V  =  V  +  V'  ,  p  =  p  +  p  '  ,  p  '  s  P 


—2  U 1 

'/a  ,  —  << 


<<  1. 


For  the  boundary  conditions  given  by  Eq.  (8),  it  is  possible  to 
show  that 

0  =  U’p"  a/P' 

and  (  5^ 

0  ?  (l+0)/(l-0). 

It  is  known  that  both  0  and  0  are  real  and  independent  of  the 
frequency  of  oscillation.  From  the  definition  of  P^-  ,  we  have 

( y  -  1 )  U2  y/y  -  1 

Pt  =  P  (1  +  2  a2  )  •  (5E 

For  small  M,  the  above  equation  can  be  approximated  by 


P  r  Pt  (  1  2RT  ^  * 


In  terms  of  the  small  adiabatic  oscillations,  Eq .  (56)  becomes 


p  *  Pt  ( 1  -  lil-  |  -  St™’  . 


Substituting  for  P' 


in  Eq .  (54) ,  we  have 


PfcUU  ’ 
RT 


p  aRT  „  _  a 

P  IJ  IT 


and 


3 


~  M  ~  1 
'  T  +  M 


1. 


(58) 


For  the  inflow  conditions  of  Eq.  (25),  0=0  and  g  =  1  and 
they  represent  adiabatic  solid  wall  boundary  conditions  involving 
specification  of  U,  V,  and  3T/9n  =  0.  The  inflow  conditions  act 

as  a  perfect  reflecting  boundary,  independent  of  M. 


The  analysis42  0f  the  inflow  condition  of  Eq .  (31)  shows 


0  =  -M  and  g  s  =  1.  (59) 

For  the  exit  boundary  conditions  (P  =  Pe )  of  Eq .  (43), 

Fedorchenko^  showed  that  Q  =  «  and  g  =  -1 ,  which  indicates  that 
the  conditions  act  like  a  perfect  reflecting  boundary. 

For  the  total  transmission  of  the  quasi-plane  normal  wave 
through  the  inlet  and  exit  boundaries  (that  is,  |  6  j  <<  1),  it  is 
required  that  ( O ) i n  f iow  =  -1  or  (0)Outflow  =  T.  If,  however, 

inflow  >  0  or  (0)Outflow  <  0'  Tt  may  lead  to  the  amplification 
of  the  longitudinal  oscillations  in  the  domain  of  interest. 

The  foregoing  analysis  does  not  take  into  account  the  radial 

oscillations  arising  from  the  inhomogeneous  boundary  conditions 

at  the  inlet  and  the  exit  (if  t  0,  then  *  0)  ,  as  well  as 

3  r  9  r 

the  refraction  effects  in  the  vortex  flow. 42 

It  is  useful  to  note  that  the  reason  for  considering 
perfectly  reflecting  permeable  boundary  conditions  at  the  exit 
(except  for  Eq.(37)  from  Rudy  and  St r i kwe rda j  is  du0  to 
Hankey,46  who  suggested  the  need  for  a  feedback  mechanism  to 
sustain  the  inherent  unsteady  behavior  which  is  characteristic  of 
the  bluff-body  shedding  phenomenon  at  high  Reynolds-number  flows. 
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3  . 


COMPUTATIONAL  DETAILS 


Using  the  time-dependent,  ax isymmetr ic ,  compressible 
Nav ier-Stokes  equations  of  a  perfect  gas,  solutions  are  obtained 
for  the  subsonic  flowfield  in  the  CBCC.  Figure  40  shows  schema¬ 
tically  the  APL  CBCC  and  the  physical  domain  which  consists  of  60 
axial  and  46  radial  nodes.  Although  the  dynamic  behavior 
observed  in  the  centerbody  near  wake^-S  involved  combustion  due  to 
the  annular  air  stream  and  central  fuel  jet,  present  computations 
are  directed  towards  obtaining  a  physically  meaningful  unsteady 
solution  involving  the  flowfield  due  to  the  annular  stream  alone. 

The  equations  describing  the  conservation  of  mass,  momentum, 
and  energy  may  be  written  in  the  following  conservative  form  in 
physical  space  as 


3_E  3  F  1  3  rG 

8t  (Tx  r  a  r 


H 


whe  re 


P 

PU 

PU 

pU2  -  T  xx 

E  = 

P  v 

,  F  = 

p  UV  -  t  rx 

P  e 

p  Ue  -  Ut  xx  -  Vxxr  -  KTX 

P  v 

0 

p  UV  -  t  rx 

0 

p  V2  -  t  rr 

,  H  = 

"  T00/r 

pVe  -  Vt  rr  -  Utxr  -  kTr 

0 

(60) 


(61) 


89 


and  where 


x  =(2u+A)tH  +  Xf-  +  ±-i)  -  P 

xx  (Tx  {  r  3  rJ 


i  =  ( 2\x  +  A  )  jLX  + 
rr  5r 


V  ,  3  U 


)  -  P 


(2m  +  A  )  —  +  A  ( 


3  U  ,  3  V  \ 

rx  yfJ 


T  =  T 

xr  rx 


3  U  3  V\ 


The  values  of  Prandtl  number  (Pr)  =  0.72  and  ratio  of  speci¬ 
fic  heats  (y  )  =  1.4  are  specified.  The  molecular  viscocity  u  is 

related  to  the  temperature  by  Sutherland's  law.  Finally,  the 
pressure  is  related  to  temperature  by  the  equation  of  state 

P  =  p  RT.  ( 66 ) 

These  equations  are  advanced  in  time  using  the  vectorized 
computer  code  of  Shang^7  which  employs  MacCormack's  explicit  and 
unsplit  algorithm. 4^  A  forward-predictor,  backward-corrector 
operator  sequence  is  used  for  each  time  step.  Exponential  grid 
str  telling  is  used  to  adequately  resolve  the  flowfield  in  the 
anticipated  regions  of  large  gradients.  The  physical  domain  is 
transformed  into  a  computational  domain  of  unit  square  with  equal 
spatial  step  sizes,  Ax  and  Ar  (Ax  *  fir)  .  The  time  step  of 
compressible  equations  is  limited  by  the  Courant-Fr iedr ichs-Lewy 
(CFI.)  condition.  The  maximum  allowable  time  step  in  transformed 
coordinates  is  given  by 


T 

A  T  =  1  /  -- 
CKL  [AA 


+  AT  + 


i^x  ,  ’  x  .  2 

a<T  +  nr>  + 


r,  2,  ,  1/2 


w  h  o  r  • 


-  r  i 
’  x 


t-  r  V 
r 


M,  ■=  ••  n  f  •:  v 

f  (x)  and  r  (r)  are  the  transformed  variables  needed  for  the 
t  r  i  ns f ormn t i on  of  the  physical  domain  to  the  computational 
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4Gb.  Temporal  Variation  in  Combustor  Mass  Flow  and  AT  ^  f°r  caso  3b  of  Table 
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Figure  43b.  Temporal  Variation  in  Combustor  Mass  Flow  and  AT  for  Case  2  of  Table 
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terbody  diameter  Dc)  are  equal  to  1.1  and  0.5  respectively.  It  is 
observed  that  the  initial  development  of  the  flowfield  is  accom¬ 
panied  by  the  generation  of  longitudinal  harmonic  oscillations 
with  a  clearly  defined  dominant  frequency  wx]^  ~  182  Hz  which  is 
close  to  the  fundamental  characteristic  frequency  of  a  one¬ 
dimensional  quarter-wave  resonator  tox^  ~  a/^Ld  191.  cycles/sec. 

In  addition,  resonance  also  occurs  at  all  harmonics  of  the  fun¬ 
damental  frequency,  i.e.,  wX2  =  2uixl'  •  ••>  wxn  =  ruo^.  The 
longitudinal  oscillations  attenuated  in  about  3,500  time  steps 
which  clearly  corresponds  to  the  time  it  takes  for  a  particle 
entering  the  inflow  boundary  to  reach  the  exit  plane, 

Ld/uinflow  ~  8  milliseconds.  The  vorticity  contours  of  Figure  42 
show  heavy  concentration  in  the  boundary  layers  of  the  duct  wall 
and  the  centerbody  wall.  The  positive  vorticity  contours 
corresponding  to  counterclockwise  vorticity  originating  from  the 
duct  wall  boundary  layer  are  shown  by  hatching.  The  negative 
vorticity  contours  originating  at  the  centerbody  boundary  layer 
are  carried  past  into  the  centerbody  near  wake  by  convection.  It 
is  observed  that  the  tail-end  of  the  zero-vorticity  contour 
extends  in  the  axial  direction  towards  the  outflow  boundary  at  a 
rate  of  3  m/s.  The  stretching  of  the  recirculation  region  is  the 
only  unsteady  feature  of  the  flowfield  after  the  combustor  mass 
flow  reaches  the  steady  state. 

Figure  43  of  Case  2  shows  that  there  is  no  variation  in  the 
mass  flow  rate  or  in  AT(-pL  °f  the  flowfield  with  time.  The 
pressure  contour  plot  seen  in  Figure  44  clearly  demonstrates  the 
existence  of  a  shear  layer  originating  from  the  trailing  edge  of 
the  centerbody  and  shows  the  location  of  the  vortex  center  in 
this  shear  layer.  Furthermore,  as  may  be  anticipated,  the 
pressure  minimum  occurs  at  the  vortex  center  and  the  pressure 
maximum  at  the  centerline  stagnation  point  of  the  centerbody. 

All  other  saliert  features  of  the  flowfield  are  identical  to  the 
previously  discussed  Case  1. 

Figures  15  to  47  correspond  to  Cases  3a  and  3b.  Case  3a 
corresponds  to  a  mass  flow  rate  of  2.5  kg/s  and  a  Reynolds  number 
of  105  on  a  (5  0  x  45)  grid  and  Case  3b  corresponds  to  a  mass  flow 
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ure  42.  Velocity-Vector  and  Vort ic ity-Contour  PIo 
1  of  Table  7. 
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Reynolds  number  based  on  the  duct  diameter  and  the  average  Mach 
number  of  the  flowfield  based  on  the  average  inflow  velocity. 
Table  7  also  provides  the  details  of  the  number  of  timesteps 
(iterations)  and  the  corresponding  total  time  in  seconds  over 
which  the  flowfield  is  computed.  Quantities  shown  in  parentheses 
in  Table  7  must  be  interpreted  with  care.  The  asterisk  denotes 
the  data  of  Reference  39. 

The  numerical  calculations  are  performed  on  the  AFWAL  CDC 
CYBER  750/175  computer  and  on  the  NASA  Lewis  Research  Center 
CRAY-1  computer.  The  original  AFWAL  computer  code  has  been 
optimized  for  the  CBCC  con f igura t ion  by  reducing  the  field-length 
of  computer  memory  from  the  initial  requirement  of  363,0000  words 
to  172,100^  words  for  the  baseline  domain.  The  optimized  ver¬ 
sion  of  the  computer  code  requires  a  field-length  of  216,2400 
words  for  the  extended  domain. 

The  relative  computational  efficiency,  or  the  so-called 
"computational  effort"  t  is  defined  as  the  CPU  time  required  to 
advance  the  solution  by  one  time  step  per  spatial  grid  point; 
i.o.,  t  =  CPU  in  seconds/number  of  mesh  points  x  number  of  time 
steps.  For  the  present  computer  code  x  -  5.96  x  10~4  seconds  on 
the  CYBER  and  r  =  4.4  x  10“5  seconds  on  the  CRAY.  All  the 
numerical  computations  are  performed  for  the  (60  x  46)  grid  on 
the  CYBER  and  for  the  (80  x  46)  grid  on  the  CRAY.  For  the  solu¬ 
tion  to  advance  10,000  time  steps,  the  approximate  CPU  time 
required  for  the  computer  code  is  4.5  hours  on  the  CYBER  and  0.45 
hours  on  the  CRAY. 

The  numerical  calculations  of  Cases  1  through  4  of  Table  7 

show  that  the  flowfield  attains  a  steady  mass-flow  rate  and  the 

time  required  for  the  flowfield  to  reach  this  converged  mass-flow 

rate  varied  anywhere  between  L^/U  j_nf  ]_ow  to  IL^/U  iow.  The 

velocity  vector  plot  of  Case  1,  Figure  42,  shows  a  stationary 
recirculation  region  in  the  centerbody  near  wake  am!  does  not 
indicate  any  evidence  of  vortex  shedding.  The  r  nterline  reat¬ 
tachment  point  extends  all  the  way  to  the  exit;  plane.  The  axial 
and  radial  coordinates  of  the  vortex  center  (normalized  by  cen- 
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domain.  The  calculations  are  performed  using  third-order  artifi¬ 
cial  damping  parameters  of  Reference  44  to  suppress  spurious 
numerical  oscillations. 

4.  RESULTS  AND  DISCUSS  TON 

Table  6  shows  the  details  of  nine  cases  modeled  with  de¬ 
ferent  combinations  of  several  inflow  ami  out  (low  bounder/  con¬ 
ditions.  To  identify  the  effect  of  the  relative  location  of  f  he 
inlet  and  exit  boundaries  on  the  reflection  character i s t i cs  of 
the  boundary  conditions  to  be  examined,  two  different  finite- 
difference  domains  are  considered.  The  baseline  domain 
corresponding  to  Figure  40  consists  of  60  nodes  in  the  axial 
direction  and  46  nodes  in  the  radial  direction.  '['he  domain 
represents  a  finite-cylindrical  duct  of  length,  Lg  =  45 .30  cm  and 
radius,  Rg  =  12.70  cm  and  a  centerbody  of  lengih,  Lc  =  10.60  cm 

and  radius,  Rr  =  7.00  cm.  The  mesh  is  constructed  with  variable 
step  sizes,  Ax  and  A r  to  provide  a  finer  mesh  at  the  near-wall 
regions  of  the  centerbody  and  the  duct,  as  well  as  in  the  recir¬ 
culation  zone  of  the  centerbody  near  wake.  Exponential, 
stretching  is  used  to  construct  the  mesh  with  Ax  varying  between 
a  maximum  of  1.4  cm  to  a  minimum  of  0.24  cm  and  Ar  varying 
between  0.605  cm  to  0.24  cm.  Figure  41  shows  the  extended 
domain,  which  consists  of  80  nodes  in  the  axial  direction  and  46 
nodes  in  the  radial  direction.  The  physical  dimensions  of  the 
two  domains  differ  only  in  the  axial  length  of  the  cylindrical 
duet  (  Lg  is  equal  to  67.06  cm  for  the  extended  domain). 

Table  7  shows  the  details  of  the  numerical  computations  for 
model ing  of  the  nine  different  cases  of  inflow  and  outflow  boun¬ 
dary  conditions  given  in  Table  6.  The  table  includes  the  values 
of  the  parameters  to  be  specified  for  the  inflow  and  outflow 
boundary  conditions,  CFL  number,  damping  pa  name  tors  (Cx  and  Cr) 
and  the  initial  distribution  of  the  flowfieio.  In  all  the 
numerical  computations,  the  values  of  the  temperature  and  com¬ 
ponents  of  velocity  to  be  specified  at  the  centerbody  and  one i 
walls  are  kept  constant  at  Tw  =  63b. 6°R  and  nw  -  Vw  =  0.  Table  7 

also  includes  the  results  of  the  numerical  c-  :■  ,  stations  in  terms 
of  the  average  mass  flow  of  the  eoiabus  tor  (  'awl  i  <.*  Id,  average 
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Figure  47b.  Vort ic i ty-Contour  Plots  for  Case  3b  of  Table 


ressure-Contour  Plots  for  Case 


rate  of  5.0  kg/s  and  a  Reynolds  number  of  2  x  10^  on  a  (80  x  46) 
grid.  Neither  the  increased  Reynolds  number  nor  the  variation  in 
the  relative  location  of  the  inflow  and  outflow  boundaries  show 
any  difference  from  the  results  of  Cases  1  and  2.  However,  there 
are  minor  fluctuations  in  the  mass  flow  rate  with  time, 
corresponding  to  a  fundamental  frequency  of  =  a/4L^  =  128  Hz 

for  the  numerical  calculations  performed  on  the  (80  x  46)  grid  of 
Case  3b. 


For  the  inflow  boundary  conditions  involving  the  specifica¬ 
tion  of  total  pressure,  total  temperature,  V  =  0,  or  =  0, 

it  becomes  practically  impossible  to  specify  at  the  inflow  cross- 
section  arbitrary  profiles  of  U  with  characteristic  values  of  the 
transverse  gradients  (especially  with  a  definite  value  of  mass 
flow  at  the  inflow).  When  these  inflow  conditions  are  used  in 
conjunction  with  the  outflow  conditions  involving  P  =  Pe ,  or  the 
nonreflecting  condition  of  Rudy  and  Strikwerda43  for  the 
flowfield  with  a  steady-state  solution,  it  becomes  extremely  dif¬ 
ficult  for  the  system  to  attain  the  desired  mass  flow  rate.  This 
difficulty  is  evident  from  the  results  of  Cases  1,  3a,  and  3b, 
where  the  average  mass  flow  rates  attained  by  the  flowfield  are 
given  by  2.3  kg/s,  2.5  kg/s,  and  5.0  kg/s,  respectively.  The 
situation  becomes  even  worse  for  unsteady  flowfields  or  when 
these  inflow  conditions  are  used  in  conjunction  with  outflow  con¬ 
ditions  involving  the  split  condition  of  Hasen,44  which  permits 
variation  of  pressure  with  time  at  the  outflow  boundary. 

The  initial  distribution  of  the  flowfield  for  Case  4 
corresponds  to  the  solution  of  Scott  and  Hankey^  at  25,700  DT. 
Figure  48  shows  that  the  fluctuations  in  the  mass  flow  persist  for 
the  first  8,000  DT.  For  the  next  14,000  DT ,  both  inflow  and 
outflow  mass  flows  reach  essentially  a  constant  value  of  5.5  kg/s 
but  with  a  slightly  increasing  monotonicity.  The  actual  dynamics 
of  this  tendency  to  steady  mass  flow  rate  can  be  seen  from  the 
velocity  vector  plots  of  Figures  40a  and  49b  and  the  vorticity 
contour  plots  of  Figures  49c  and  49d.  In  Figures  49c  and  49d, 
breaking  up  of  the  clockwise  recirculation  vortex  occurs  as  a 
counterclockwise  vortex  originates  at  the  centerline  of  the  cen- 
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4 8b.  Temporal  Variation 


5. SO 


112 


iation  in  Combustor  Mass  Flow  for  Case  4  of  Table 


Figure  49a.  Velocity  Vector  Plots  for  Case  4  of  Table  7. 
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N  =  22,000 
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Figure  49b.  Velocity  Vector  Plots  for  Case  4  of  Table  7. 
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N  =  22,000 


Figure  49d.  Vort ic ity-Contour  Plots  for  Case  4  of  Table  7. 
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terbody  near  wake,  gains  strength,  and  grows.  At  5,000  DT,  the 
figures  show  the  remnant  of  the  clockwise  vortex  which  is  about 
to  be  washed  off  downstream.  From  this  point  onwards  the  mass 
flow  rates  tend  to  their  constant  values,  the  clockwise  recir¬ 
culation  vortex  begins  to  grow,  pushing  the  eddies  towards  the 
exit,  and  the  final  state  resembles  that  of  Figures  49b  and  49d 
depicting  the  stretching  of  a  stationary  recirculation  vortex, 
extending  in  the  axial  direction  at  a  rate  of  3m/sec.  The  large 
mass  flow  fluctuations  of  Figure  48  correspond  to  the  time  during 
which  eddies  are  washed  off  the  exit.  Thus,  approximately  8,000 
DT  are  required  for  the  initial  transient  fluctuations  to  disap¬ 
pear  and  the  flow  to  attain  a  steady  mass  flow  rate  of  5.5  kg/s. 

It  is  important  to  note  that  none  of  the  calculations 
involving  Cases  1  to  4  of  Table  7  are  continued  long  enough  for 
the  reverse  flow  to  reach  the  outflow  boundary,  which  would  even¬ 
tually  happen  due  to  the  stretching  of  the  recirculation  zone  in 
the  axial  direction.  The  rate  at  which  the  recirculation  zone 
extends  in  the  axial  direction  is  found  to  be  very  slow  at  3  m/s 
in  all  the  cases  considered. 

Figures  50  and  51  correspond  to  Case  5  of  Table  7.  The  in¬ 
flow  and  outflow  conditions  are  identical  to  Case  4,  except  for 
the  inflow  boundary  conditions,  which  are  programmed  by  the 
simplified  procedure  described  in  Paragraph  III.2.a(2b).  During 
the  first  20,000  DT,  after  the  initial  attenuation  of  the  longi¬ 
tudinal  oscillations  a  steady  mass  flow  rate  is  established  as 
evidenced  by  Figures  5la,  51c,  and  51e.  Recall  that  three  boun¬ 
dary  conditions  are  to  be  specified  when  there  is  a  reverse  flow 
at  the  outflow  boundary.  As  the  extending  recirculation  vortex 
causes  reverse  flow  at  the  outflow  boundary  around  22,000  DT , 
the  split  boundary  condition  involving  the  specification  of  only 
exit  pressure  cannot  objectively  handle  the  appropriate  physics 
of  the  flowfield.  As  soon  as  this  happens,  reverse  longitudinal 
oscillations  are  established,  causing  the  domain  to  act  as  a 
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Temporal  Variation  in  Combust or  Mass  Flow  and  AT  for  Case  5  of  Table 


Figure  50b.  Temporal  Variation  :'n  Combustor  Mass  Flow  and  ATCFL  for  Case  5  of  Table 
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Velocity-Vector  Plots  for  Case  6b  of  Table 
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resonator  with  a  fundamental  frequency  of  the  quarter-wave  reso¬ 
nator.  This  process  of  acoustic  interaction  with  the  recir¬ 
culation  zone  is  accompanied  by  shedding-like  phenomena  seen  in 
Figures  51b,  51d,  and  5 1 f .  The  mass-flow  fluctuations  of  Figure 

50  show  that  the  amplification  of  the  longitudinal  oscillations 
result  in  a  complete  reversal  of  the  mass  flow  at  the  inflow  and 
outflow  boundaries  and  the  numerical  solution  becomes  unstable 
and  nonphysical. 

Cases  6a,  6b,  and  6c  of  Table  7  refer  to  the  inflow  and 

outflow  conditions  discussed  in  Paragraphs  III.2.a(2b) 
and  III.2.b(3b).  Figures  52  and  53  are  taken  from  Reference  39. 
Figure  52  shows  that  the  mass  flow  varies  between  a  maximum  of 
4.875  kg/s  to  a  minimum  of  3.125  kg/s  and  the  variations  in  the 
ATcfl  a  result  of  the  temperature  fluctuations  in  the 
flowfield  as  ATgpp  ~  1//T*.  The  numerical  calculations  of  the 

isothermal  flowfield  showed  temperature  fluctuations  ranging  from 
4 , 000  °  R  to  250  °  R  but  the  fluctuations  are  restricted  to  a  very 
narrow  band  of  the  grid  near  the  centerline  of  the  reattachment 
point  of  the  recirculation  region.  It  appeared  that  these  tem¬ 
perature  fluctuations  were  caused  by  a  violation  of  the  CFL  cri¬ 
terion  due  to  an  earlier  discrepancy  in  the  TIMESTEP  subroutine. 
Figure  53  shows  the  velocity  vector  plots  and  vorticity  contour 
plots.  For  a  discussion  of  the  results  on  the  unsteady  solution 
of  Case  6a,  the  reader  is  referred  to  Reference  39. 

For  Case  6b,  the  numerical  calculations  are  performed  on  a 
(80  x  46)  grid  with  the  corrected  TIMESTEP  subroutine.  Figure  54 
shows  that  the  mass  fluctuations  persist  for  as  long  as  36,000  DT 
befor  establishing  a  steady  mass-flow  rate.  However,  the  tem¬ 
poral  variation  of  AT^p^  and  the  associated  temperature  fluc¬ 
tuations  (AT^pp  ~  1  //~T~)  are  indications  of  nonphysical  behavior 

of  the  solution  as  the  calculations  are  performed  under  isother¬ 
mal  conditions.  Figures  55a  to  55e  show  the  initial  oscillatory 
flowfield  before  the  establishment  of  a  single  recirculation 
vortex . 
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Pressure-Contour  Plots  for  Case  5  of  Table 


Contour  Plots  for  Case  5  of  Table 


Figure  51b.  Velocity-Vector  Plots  for  Case  5  of  Table 
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Figure  55c.  Vort ic ity-Contour  Plots  for  Case  6b  of  Table 


Figure  55d.  Vort icity-Contour  Plots  for  Case  6b  of  Table 


Figure  55e.  Pressure-Contour  _oi. 


Numerical  calculations  o£  Case  6c  are  performed  at  a  CFL 
number  of  0.8,  with  all  other  conditions  being  identical  to  Case 
6b.  Figure  56  shows  that  the  mass  flow  rate  remains  stable 
around  4.75  kg/s  except  for  the  fluctuations  resulting  from  the 
eddies  being  washed  off  the  outflow  boundary  starting  at  28,000 
DT  and  38,800  DT.  However,  the  large  temperature  fluctuations 
observed  (AT^p^  ~  1//T)  result  from  the  particular  choice  of 
boundary  conditions.  Figure  57  shows  the  velocity  vector  plots 
and  vorticity  contour  plots.  The  numerical  solutions  of  Cases  6a 
to  6c  show  that  the  character  of  mass  flow  fluctuations  and  of 
the  associated  oscillatory  behavior  of  the  flowfield  is  a  func¬ 
tion  of  CFL  number,  as  well  as  the  relative  location  of  the 
inflow  and  outflow  boundaries.  More  important,  the  solution  is 
nonphysical  because  of  the  large  temperature  fluctuations  asso¬ 
ciated  with  the  flowfield. 

Cases  7a  and  7b  are  modeled  with  the  inflow  conditions  of 
Paragraph  III.2.a(4)  and  outflow  conditions  involving  split  con¬ 
dition  of  Hasen44  [see  Paragraph  II I . 2 . b ( 3a ) ) .  However,  the 
split  condition  is  modified  for  Case  7b  as  seen  below: 

P"^1 .  =  P  if  j  <  43  and  =  0  if  j  >  43.  (68) 

The  condition  3  P/9 x  =  0  is  restricted  to  the  duct-wall  boundary 
layer.  The  advantage  of  these  boundary  conditions  is  that  they 
allow  for  the  specification  of  inflow  mass  flow  and  have  the 
desired  reflection  characteristics  for  the  modeling  of  the 
unsteady  flow  behavior. 37 

Figures  58  and  59  show  that  the  appearance  of  the  reverse 
flow  at  the  outflow  boundary  triggers  mass  flow  oscillations 
starting  with  a  fundamental  frequency  of  ~  128  Hz.  The 

establishment  of  the  longitudinal  oscillations  is  preceded  by 
fluctuations  in  AT^FL  for  a  short  duration,  starting  from  17,700 
DT  to  19,700  DT.  This  appears  to  resemble  the  beating  phenomenon 
of  two  harmonic  waves  of  different  amplitude  and  slightly  dif¬ 
ferent  frequencies.  The  thermal  waves  are  restricted  to  a  narrow 
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Figure  57c.  Vort ic ity-Contour  Plots  for  Case  6c  of  Table 


Figure  57e.  Pressure-Contour  Plots  Eor  Case  6c  of  Table 
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Figure  58a.  Temporal  Variation  o£  Combustor  Mass  Flow  and  AT  ^  £°r  caSG  Ta^e 
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Figure  62c.  Temporal  Variation  of  Combustor  Mass  Flow  and  A ^CFL 
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Figure  62b.  Temporal  Variation  oE  Combustor  Mass  Flow  and  AT  for  Case  8a  o£  Table 
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igure  61c.  Vorticity-Contour  and  Velocity-Vector  Plots  for  Case  7b  of  Table 


Figure  60b.  Temporal  Variation  in  Combustor  Mass  Flow  and  AT  for  Case  7b  of  Table 
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band  of  grid  nodes  near  the  reattachment  point  of  the  recir¬ 
culation  region.  Figure  60  of  Case  7b  shows  that  the  flow 
appears  to  be  steady  for  the  first  50,000  DT  followed  by  tem¬ 
perature  fluctuations  even  before  the  mass  flow  variations  begin. 
Figure  61  shows  the  extending  recirculation  region  before  the 
numerical  solution  becomes  unstable  and  nonphysical  due  to  the 
reverse  flow  at  the  outflow  boundary. 

The  outflow  boundary  conditions  of  Cases  8a  and  8b  involve 
specification  of  three  variables,  when  there  is  a  reverse  flow  at 
the  boundary  (see  Paragraph  IIT.2.b(5).  Unlike  the  numerical 
solutions  of  Cases  5  to  7  we  studied,  the  AT^pL  of  Figure  62 
remains  invariant  with  time  as  it  should  be,  even  if  the 
fiowfield  becomes  unsteady.  When  the  reverse  flow  appears  at  the 
outflow  boundary  around  54,000  DT ,  longitudinal  oscillations  are 
established  at  a  definite  fundamental  frequency  of  =  128  Hz 

with  a  time  period  of  7.02  milliseconds.  This  corresponds  to  the 
quarter-wave  frequency  of  the  resonator.  As  the  time  progresses, 
it  excites  the  adjacent  harmonics  of  frequency  2uix^  . .., 

nu)xi.  The  bimodal  frequency  of  the  mass  flow  oscillations 
starting  from  62,000  DT  clearly  appears  to  be  qualitatively  simi¬ 
lar  to  the  wave  shown  in  Figure  63,  which  is  a  superposition  of 

two  simple  harmonic  waves,  {  (2.0  -  sin  ut  +  (0.5  -  ^-)  sin  3wt} 

The  two  dominant  frequencies  of  the  fiowfield  correspond  to 
w x ]i  and  o)3xl*  However,  it  takes  approximately  seven  cycles 
for  the  longitudinal  oscillations  to  attenuate  before 
establishing  a  steady  mass-flow  rate.  The  velocity  vector  plots 
and  vorticity  contour  plots  for  the  first  40,000  DT  of  Figures  61 
and  64  are  identical  and  are  not  repeated  in  Figure  64.  Figures 
65  and  66  correspond  to  inflow  mass  rate  of  4  kg/s  of  Case  8b 
and  the  results  are  similar  to  Case  8a. 

Figure  67  shows  the  velocity-vector  and  vor t ic i ty-contour 
plots  of  Case  9.  The  outflow  boundary  conditions  of  Reference  45 
allows  for  the  reverse  flow  to  appear  at  the  outflow  boundary 
without  causing  any  oscillations  in  the  mass  flow  rate  but  the 
results  do  indicate  high  temperature  fluctuations. 
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59c.  Vort icity-Contour  Plots  for  Case  7a  of  Table 


Figure  59b.  Velocity-Vector  Plots  for  Case  7a  of  Table 


59a.  Velocity-Vector  Plots  for  Cas 
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Figure  64a.  Velocity-Vector  Plots  for  Case  8a  of  Table 
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Figure  66d.  Vort ic ity-Contour  Plots  for  Case  8b  of  Table 


Figure  67.  Velocity-Vector  and  Vor t ic i ty-Contour  Plo 
9  of  Table  7. 
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SECTION  IV 

CONCLUSIONS  AND  RECOMMENDATIONS 


This  section  presents  the  main  conclusions  emerging  from  the 
numerical  flowfield  modeling  of  the  centerbody  combustor 
configuration,  followed  by  recommendations  for  further  activity. 

1.  CONCLUSIONS 

Numerical  predictions  of  the  turbulent  recirculating  flows 
in  the  centerbody  combustor  configuration  have  been  made  for 
nonreacting  situations  within  the  frameworks  of  both  the 
Reynolds-averaged  and  time-dependent  Nav ier-Stokes  equations. 

The  steady-state  predictions  with  the  Reynolds-averaged  equations 
have  examined  both  the  large-scale  and  the  small-scale  CBCC ,  as 
well  as  the  confined  two-dimensional  mixing  layer.  The  turbulent 
eddy  viscosity  is  computed  by  the  k-e  model  and  the  calculations 
considered  the  effect  of  streamline  curvature  and  the 
preferential  influence  of  the  normal  stresses  on  the  eddy 
viscosity.  For  the  CBCC  flowfields,  the  computational 
investigations  address  both  the  annular  (air)  and  central  (CO2) 
flows.  The  predicted  results  of  the  axial  and  radial 
distributions  of  the  mean  and  fluctuating  axial  and  radial 
velocity  components  are  compared  with  the  available  experimental 
data . 


The  time-dependent  calculations  of  the  CBCC  address  the 
large-scale  configuration  with  only  the  annular  (air)  stream 
present.  No  modeling  of  turbulence  is  included  in  the  present 
calculations  and  the  formulation  accounts  for  the  molecular 
viscosity  only  (except  for  the  effects  of  artificial  viscosity  in 
the  numerical  computations).  The  predictive  calculations  with 
the  time-dependent  equations  have  examined  the  influence  of  the 
boundary  conditions  at  the  inlet  and  exit  on  the  interior 
flowfield  from  the  viewpoint  cf  explaining  the  presence  or  the 
absence  of  a  dynamic  behavior. 
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Following  are  the  major  conclusions  of  this  study: 

n  The  numerical  calculations  of  the  Reynolds-averaged 

Nav ie r-Stokes  equations  with  the  k-c  turbulence  model 
provide  qualitatively  correct  predictions  that  are 
consistent  with  experimental  observations  in  the  confined 
two-dimensional  mixing  layer  and  the  CBCC  flowfields. 

»  The  quantitative  predictions  in  these  configurations  are 
generally  good  for  the  mean  field  and  poor  for  the 
fluctuating  field.  This  observation  is  characteristic  of 
many  complex  turbulent  flows  and  is  largely  due  to  the 
inadequacy  of  the  standard  k-e  model  (since  it  is  an 
isotropic  model  and  also  does  not  account  for  the  extra 
strain  rates  arising  from  interacting  shear  layers,  large 
streamline  curvature,  etc.). 

•  Ad  hoc  modifications  to  the  standard  k-t:  model  to  account 
for  the  effects  of  large  streamline  curvature  and  the 
preferential  influence  of  normal  stresses  emphasize  the 
complex  nature  of  the  CBCC  flowfields  by  exhibiting  a 
parametric  dependence  on  the  geometric  scale  and  flow  rates. 
For  example,  the  predictions  of  the  small-scale  CBCC  with 
the  standard  k-c  model  and  the  improvements  thereon  due  to 
streamline-curvature  correction  and  dissipation-equation 
modification  conform  well  to  the  behavior  observed  in 
several  recirculating  flows  by  other  investigators.  The 
predictions  of  the  large-scale  CBCC,  however,  exhibit  trends 
that  are  considerably  different,  since  the  standard  model 
appears  to  be  the  best  overall,  while  the  modifications 
overestimate  the  recirculation-zone  lengths. 

•  The  predictions  of  the  vortex-center  characteristics  in  the 
small-scale  and  large-scale  CBCC's,  on  the  other  hand, 
demonstrate  similar  trends  vis-a-vis  the  experimental  data. 
Since  the  underprediction  of  the  axial  location  of  the 
vortex  center  by  the  standard  model  is  consistent  with  the 
experience  of  other  investigators,  it  suggests  that  either 
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the  nature  of  the  prediction  is  different  in  different 
regions  of  the  flowfield,  or  the  measurement  of  the  cen¬ 
terline  recirculation  length  in  the  large-scale  CBCC  is 
significantly  affected  by  the  LDA  technique. 

•  Preliminary  results  appeared  to  demonstrate  the  capability 
of  the  computations  with  time-dependent  Nav ier-Stokes 
equations  to  simulate  the  unsteady  flow  features  in  the 
CBCC.  The  physical  realism  of  the  predicted  behavior  is  not 
clear,  however,  in  view  of  the  large  axial  extent  of  the 
recirculation  region  and  the  very  large  fluctuations  in  the 
mass  flow  rates  and  temperature.  One  possible  source  of 
this  difficulty  is  the  lack  of  an  adequate  turbulence  model 
in  the  time-dependent  computations. 

•  Comprehensive  examination  of  the  inflow  and  outflow  boundary 
conditions  in  the  time-dependent  computations  suggests  an 
apparent  lack  of  uniqueness  in  the  numerical  solutions  of 
the  unsteady  flow.  Several  computational  case  studies 
appear  to  show  that  the  time-dependent  formulation  yields  a 
steady-state  solution  asymptotically.  The  fluctuations  in 
the  mass  flow  persist  for  a  few  thousand  time  steps  before 
the  mass  flow  becomes  steady.  It  remains  unchanged  for  some 
time  and  the  fluctuations  start  to  appear  again.  The  only 
unsteady  feature  in  the  flowfield  during  the  time  the 
mass-flow  rate  remains  steady  is  the  slow  axial  propagation 
of  the  recirculation  region  toward  the  exit  boundary  (at  a 
speed  of  approximately  10  ft/s).  This  process  continues 
until  the  reverse  flow  reaches  the  exit  boundary,  thereby 
rendering  the  exit-boundary  conditions  incompatible.  This 
inadequacy  of  the  outflow  boundary  conditions  to  handle  the 
reverse  flow  causes  the  flowfield  to  become  unsteady  and 
nonphysical,  breaks  up  the  rec i rcu la t ion  zone,  and  exhibits 
a  shedd i ng- 1 i ke  phenomenon  as  longitudinal  oscillations  in 
the  mass-flow  rate  begin  to  appear  with  a  definite 
fundamental  frequency  corresponding  to  a  quarter-wave 
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resonator.  It  takes  a  considerable  time  for  these 
longitudinal  oscillations  to  attenuate  before  the  flowfield 
attains  a  steady  state.  The  whole  process  repeats  again 
with  the  axial  migration  of  the  recirculation  region  as  the 
calculations  proceed  further.  Thus,  it  appears  that  the 
shedding-1 ike  behavior  observed  during  the  time  period 
associated  with  the  mass-flow  fluctuations  may  be  neither 
self-sustaining,  nor  physical.  Therefore,  a  conclusive 
demonstration  of  the  successful  simulation  of  unsteady 
features  of  subsonic  internal  flows  characteristic  of  the 
CBCC  is  not  available  at  present. 

2.  RECOMMENDATIONS 

In  view  of  the  foregoing  conclusions,  the  present  study 
offers  the  following  recommendations  for  further  activity  in  the 
time-averaged  and  time-dependent  computations  of  turbulent  mixing 
in  subsonic  internal  flowfields. 

•  For  predicting  complex  turbulent  flows  wherein  large-scale 
unsteadiness  is  not  significant,  the  time-averaged 
formulation  (involving  the  Reynolds-averaged  equations  in 
constant-density  flows  and  the  Favre-averaged  equations  in 
variable-density  flows,  especially  reacting  flows)  should  be 
adequate,  provided  higher-level  turbulence  models  are  used. 
For  flowfields  wherein  anisotropy  is  not  significant,  the 
standard  k- e  model  should  suffice.  In  the  near-wake  region 
of  the  CBCC  where  anisotropic  effects  are  important, 
turbulence  models  involving  the  direct  solution  of  the 
equations  for  the  Reynolds  stresses  are  worthwhile.  A 
useful  start  in  this  direction  would  be  the  application  of 
algebraic  stress  models. 

•  For  the  prediction  of  turbulent  flows  possessing  unsteady 
features,  the  Reynolds-averaged  formulation  which  retains 
the  explicit  time  dependence  should  be  considered.  This 
approach  should  be  adequate  for  the  numerical  simulation  of 
those  unsteady  flows  whose  characteristic  frequency  of  the 


unsteady  mean  flow  is  much  smaller  than  the  characteristic 
frequency  of  the  typical  large  eddies.  Although  it  is 
instructive  to  regard  the  present  time-dependent 
calculations  with  the  full  Navier-Stokes  equations  employing 
the  laminar  viscosity  as  a  special  subcase  of  the 
Reynolds-averaged  simulation  with  zero  turbulence,  its 
relevance  to  the  high-Reynolds-number  flows  in  the  CBCC  is 
questionable.  If  the  CBCC  flowfield  is  representative  of  a 
cyclic  unsteadiness  with  a  single  narrow-band  frequency 
(which  is  much  smaller  than  the  typical  eddy  frequency)  and 
no  additional  complications  such  as  three-dimensional  effect 
and  free-stream  turbulence  are  present,  the  axisymmetric 
Reynolds-averaged  formulation  should  prove  to  be  useful.  Of 
course,  the  earlier  remarks  regarding  the  turbulence 
modeling  still  apply. 

•  It  appears  from  the  unpublished  AFWAL/POSF  data  on  the  CBCC 
power  spectra  that  this  configuration  involves  a  broad-band 
distribution  of  frequencies.  Furthermore,  the  present 
experimental  setup  may  not  be  entirely  free  from  free-stream 
turbulence  and  three-dimensional  effect.  Thus,  it  may  be 
necessary  to  consider  the  three-dimensional  simulation  with 
the  Reynolds-averaged  equations. 

•  Nevertheless,  further  research  is  clearly  needed  to 
ascertain  the  influence  of  the  inflow  and  outflow  boundary 
conditions  on  the  successful  simulation  of  unsteady 
features,  or  the  lack  thereof  in  subsonic  internal 

f lowf ields. 

•  The  numerical  computation  of  turbulent  flows  through 
time-averaged  Navier-Stokes  equations  presumes  that  the 
instantaneous  flow  variable  can  be  decomposed  into  a 
well-defined  mean  component  (which  may  be  time-dependent) 
and  a  fluctuating  component  that  is  essentially  random. 
Presently  a  viable  theoretical  framework  to  describe 
turbulent  flows  involving  distinct  quasi-deterministic, 
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large-scale  structures  (which  occur  in  transitional  flows 
and  in  flows  which  are  younger  and  closer  to  their  origin) 
does  not  exist.  However,  the  approach  based  on  large-eddy 
simulation  with  subgrid-scale  turbulence  modeling  is 
computationally  feasible  and  theoretically  appealing,  since 
the  time-dependent,  three-dimensional  Nav ier-Stokes 
equations  are  used  to  calculate  the  eddies  larger  than  a 
selected  cutoff  length  and  turbulence  modeling  is  required 
only  for  the  eddies  that  are  smaller.  While  the  available 
computer  resources  are  not  sufficient  for  this  approach  to 
simulate  reacting  flowfields  in  the  CBCC ,  a  start  should  be 
made  with  this  approach  for  simpler  flows  which  may  be 
regarded  as  subsets  of  the  CBCC. 
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